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Abstract 


“’An  equipment  model  has  been  developed  for  the  low  pressure  chemical  vapor  deposition 
(LPCVD)  of  polycrystalline  silicon  in  a  horizontal  tube  furnace.  The  model  predicts  the 
wafer-to-wafer  deposition  rate  down  the  length  of  the  tube.  Inputs  to  the  model  include: 
silane  flow  rates  from  three  injectors,  injector  locations,  locations  of  and  temperatures  at 
three  thermocouples,  operating  pressure,  the  number  of  wafers,  wafer  diameter,  the 
location  of  the  wafer  load,  and  other  physical  dimensions  of  the  furnace  such  as  tube  length 
and  inner  diameter.  The  model  is  intended  to  aid  the  process  engineer  in  the  operation  of 
equipment,  including  the  selection  of  optimum  process  parameters  and  process  control 
based  on  measured  deposition  thicknesses.  'Fhe  model  is  also  flexible  enough  to  aid  in  the 
design  of  new  equipment. 

The  one  dimensional  finite  difference  model  encompasses  the  convective  and  diffusive 
fluxes  of  silane  and  hydrogen  in  the  annular  space  between  the  wafer  load  and  tube  walls. 
The  reaction  of  silane  is  modeled  with  full  account  taken  of  the  generation  and  transport  of 
hydrogen.  Kinetic  and  injection  parameters  in  the  model  were  calibrated  using  a  series  of 
■  nine  statistically  designed  experiments  which  varied  four  parameters  over  three  levels.  'Fhe 
model  accurately  predicts  the  axial  deposition  profile  over  the  full  range  of 
experimentation  and  demonstrates  good  extrapolation  beyond  the  range  of  experimental 
calibration.  The  model  was  used  to  predict  a  set  of  process  parameters  tiiai  wuulu  icSiili  iit 
the  least  variation  of  deposition  rate  down  the  tube.  JFie  predicted  parameters  agr'^2  well 
with  experimentally  determined  optimum  conditions.  Y 
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ABSTRACT 

An  equipment  model  has  been  developed  for  the  low  pressure  chemical  vapor  depo¬ 
sition  (LPCVD)  of  polycrystalline  silicon  in  a  horizontal  tube  furnace.  The  model 
predicts  the  wafer-to-wafer  deposition  rate  down  the  length  of  the  tube.  Inputs  to 
the  model  include:  silane  flow  rates  from  three  injectors,  injector  locations,  locai- 
tions  of  and  temperatures  at  three  thermocouples,  operating  pressure,  the  number 
of  waders,  wafer  diameter,  the  location  of  the  wafer  load,  and  other  physical  di¬ 
mensions  of  the  furnace  such  as  tube  length,  and  in:'  aiameter.  The  model  is 
intended  to  aid  the  process  engineer  in  the  operation  of  equipment,  including  the 
selection  of  optimum  process  paraumeters  and  process  control  based  on  meastired 
depositlc..*'  thicknesses.  The  model  is  also  flexible  enough  to  aid  in  the  design  of 
new  equipment. 

The  one  dimensional  finite  difference  model  encompasses  the  convective  and  dif¬ 
fusive  fluxes  of  silane  and  hydrogen  in  the  annular  space  between  the  wader  load  and 
tube  walls.  The  reaction  of  silane  is  modeled,  with  full  account  taken  of  the  genera¬ 
tion  aind  transport  of  hydrogen.  Kinetic  and  injection  parameters  in  the  model  were 
calibrated  using  a  series  of  nine  statistically  designed  experiments  which  varied  four 
paraimeters  over  three  levek.  The  model  accurately  predicts  the  axial  deposition 
profile  over  the  full  range  of  experimentation  and  demonstrates  good  extrapolation 
beyond  the  range  of  experimental  calibration.  The  model  was  used  to  predict  a 
set  of  process  parameters  that  would  result  in  the  least  variation  of  deposition  rate 
down  the  tube.  The  predicted  parameters  agree  well  with  experimentally  deter- 
mined  optimum  conditions. 
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Figure  1:  Schematic  of  LPCVD  reactor. 

1  INTRODUCTION 

1.1  Motivation  and  Background 

Low  pressure  chemical  vapor  deposition  (LPCVD)  of  polysilicon  in  tubular  hot-wall 
reactors  has  become  a  wide  spread  process  in  the  manufacture  of  semiconductor  de¬ 
vices.  The  proliferations  of  this  process  over  previous  atmospheric  cold-wall  pressure 
CVD  is  due  to  the  ability  to  process  larger  quantities  of  wafers  with  better  unifor¬ 
mities  across  wafers  and  wafer-to-wafer  uniformities  without  a  significant  loss  (5  to 
10  times  less)  in  growth  rates.  The  low  pressures  (0.2-1  torr)  increase  the  diffusion 
coefficient  by  three  orders  of  a  magnitude  making  the  process  reaction  rate  lim¬ 
ited.  A  potential  problem  is  the  accumulation  of  silane  at  the  reactor  walls  causing 
prticulates  and  undue  stress  on  the  process  tube.  This  problem  can  be  avoided  by 

adding  a  liner  inside  the  process  tube  which  can  be  easily  removed  and  replaced  for 
periodic  cleaning. 

The  BTU  Engineering/Bruce  System  7351 C  horizontal,  hot-wall  furnace  is  a 
typical  commercial  reactor  used  for  the  LPCVD  of  polysilicon,  schematically  shown 
in  Figure  1.  The  process  area  consists  of  a  quartz  process  tube  surrounded  by  a 
three  zone  heating  coil  with  a  quartz  liner  inside  the  process  tube.  The  wafer  load 
is  inserted  and  removed  on  a  silicon  carbide  cantilever  attached  to  the  front  door 
of  the  reactor.  The  wafer  load  rests  in  quartz  boats  and  is  situtated  concentrically 
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in  the  tube  so  that  the  wafers  are  perpendiculaj  to  the  main  gas  Bow.  The  vacuum 
is  pulled  from  the  source*  end  of  the  tube.  Gas  (pure  silane,  SiH4)  can  be  injected  { 
through  three  injectors  -  fixed  load  end  -  and  two  moveable  injectors  -  which  are 
usually  situated  in  the  center  and  source  cones.  The  temperature  of  the  zoned 
heater  is  controlled  during  deposition  by  three  thermocouples  which  are  situated 
in  a  quf.jtz  sheath  inside  the  process  tube  and  are  individually  located  near  the 
beginning  of  the  wafer  load,  in  the  center,  and  near  the  end  of  the  wafer  load. 

These  paraimeters  are  necessary  inputs  to  an  equipment  model. 

At  present,  processes  are  designed  by  empirical  rules.  Current  process  specificar 
tions  require  across  wafer  uniformities  on  the  order  of  ±2%  and  across  wafer  load 
uniformities  of  :^5%.  Equipment  specifications  claim  that  across  wafer  uniformities 
of  ±.1%  and  across  wafer  load  uniformities  of  ±3%  can  be  produced.  Obtaining  the 
equipment  paramter  settings  to  achieve  these  process  specifications  is  an  expensive 
and  time  consuming  chore  based  on  trial  and  error  experimentation  inside  a  small 
operating  regime  iti  which  the  equipnicjiit  ailu  process  is  well  behaved.  Investigation 
of  operating  regimes  outside  the  empirical  “envelope"  is  a  prc>hibitively  expensive 
task.  As  device  specifications  become  ever  more  stringent,  process  specifications  will 
be  more  difficult  to  meet.  Equipment  models  are  the  most  cost  effective  method  to 
achieve  the  necessary  goals  of  equipment  operators  and  equipment  designers. 

"I  O 

Many  studies  have  been  conducted  to  advance  the  physical  understanding  of  the 
LPCVD  of  polysilicon.  Hitchman  [l]  derived  a  basic  linear  model  of  the  LPCVD 
reaction  kinetics.  Van  Den  Brekel  [2]  and  Claassen  [3]  investigated  the  LPCVD  re¬ 
action  chemistry  for  polysilicon  and  developed  an  understanding  for  the  main  effects 

^Tliroaghoat  tki*  paper  the  furnace  will  be  referred  to  in  three  lonea:  the  load  sone,  front  area 
of  the  furnace  where  the  wafers  are  inserted  and  removed;  the  source  sone,  which  is  the  area  at  the 

r*%r  r^f  ihtft  unr)  thit  r#tnt«r  cnni^. 
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of  the  reacton  chemistry  as  supported  by  experimental  observations.  Middleman 
[4j  developed  a  numerical  model  for  the  mass  transport  in  an  iinnular  LPCVD  reac¬ 
tor.  He  showed  the  effect  of  including  diffusion  along  with  the  gas  convection  and 
that  the  flow  in  the  annulax  region  was  insufficient  to  create  significamt  circulating 
flows  between  the  wafers.  The  above  studies  were  invaluable  to  the  the  physical 
understanding  of  the  LPCVD  process. 

The  most  complete  equipment  models  developed  have  been  based  on  two  very 
different  premises:  numerical  solutions  of  the  system  physical  representation;  and 
expert  systems  mainly  based  on  empirical  rules.  A  numerical  model  for  the  LPCVD 
of  polysilicon  was  developed  by  Jensen  [5,6].  This  model  embodies  the  most  ad¬ 
vanced  reaction  kinetic  model  derived  from  the  studies  of  Vam  Den  Brekel  [2]  and 
Claassen  [3]  with  a  gas  flow  model  dependent  on  the  system  geometry.  The  model 
predicts  the  axial  deposition  profile  for  polysilcon  deposition  for  a  ramped  temper¬ 
ature  processes  in  which  there  is  no  gas  injection  other  than  at  the  load  end  of 
the  furnace.  It  is  important  to  note  that,  although  the  model  is  able  to  predict 
the  axial  profiles  for  which  it  was  tested,  the  model  and  data  were  not  run  with 
the  same  equipment  parameter  settings.  The  expert  sytem  approach  most  recently 
developed  at  Berkeley  [7]  for  the  LPCVD  of  doped  polysilicon.  The  expert  system 
partitions  the  deposition  process  goals  into  six  modules  for  determination  of  the 
resitivity,  thickness,  uniformity,  grain  size,  film  stress,  and  a  support  module.  £ach 
individually  search  a  data  base  of  empirical  rules  for  the  correct  equipment  settings 
to  meet  the  process  specifications.  The  expert  system  is  limited  to  operation  in 
a  small  operating  window  in  which  it  has  empirical  knowledge  and  would  benefit 
from  an  inclusion  of  equipment  noise  parameters  such  as  temperature  fluctuations 
or  gas  flow  variations  for  uniformity  analysis. 
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Figure  2:  Generic  representation  of  an  equipment  mode!. 

2  NEED  FOR  EQUIPMENT  MODELS 

2.1  Dennition  and  Configuration 

In  the  most  general  sense,  an  equipment  model  is  a  body  of  knowledge  which  pro¬ 
vides  predictions  about  the  outputs  from  a  unit  manufacturing  process,  given  in¬ 
formation  about  the  inputs  to  the  process.  Figure  2  illustrates  a  generic  equipment 
model  with  outputs  and  two  classes  of  inputs,  the  process  parameters  and  the  dis¬ 
turbances.  The  process  parameters  are  those  parioneters  that  we  exercise  direct 
control  over,  for  example,  temperature,  pressure,  gas  flow  rates.  The  dist  irbances 
are  those  inputs  to  the  process  which  are  subject  to  imintended  and  tmdesired  vari¬ 
ations.  In  some  cases,  the  magnitude  of  the  disturbances  can  be  monitored,  while  in 
other  cases  they  cannot.  Examples  of  disttirbances  include  variations  in  the  proper¬ 
ties  of  incoming  material  and  vwiations  in  the  process  parameters  themselves.  The 
alternate  labels  presented  in  parentheses  in  Figure  2,  (control  factors,  noise  factors, 
response  factors)  are  the  corresponding  terminology  tised  by  Taguchi  [8]  to  describe 
manufacturing  processes. 

At  a  minimum,  a  competent  equipment  model  must  provide  information  about 
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the  ouputs  given  information  about  the  process  paraunetcrs.  A  more  useful  model, 
will  also  include  predictions  about  the  effect  of  disturbances.  It  is  vital  that  the 
model  provide  information  about  the  variation  of  the  outputs  and  the  dependence  of 
this  variation  on  the  process  paramters  and  the  disturbances.  The  variations  might 
include  across  the  wafer  variation,  wafer  to  wafer  variation,  and  batch  to  batch 
variation.  A  competent  model  might  address  all  three  classes  cf  variation,  or  might 
focus  on  the  most  important  class  as  indicated  by  experience.  Accurate  predictions 
about  the  process  mean  are  in  general  less  important  than  accurate  predictions 
about  variation,  since  in  most  processes  the  mean  can  be  adjusted  tu  its  target 
value  without  a  substantial  defect  on  variation.  An  example  of  such  an  adjustment 
would  be  the  length  of  time  in  an  LPVCD  deposition.  Occasional  exceptions  to  this 
generalization  will  require  that  the  model  accxirately  predict  the  mean  as  well. 

2.2  Uses  for  Equipment  Models 

Equipment  models  can  be  used  to  aid  in  the  operation  of  existing  equipment,  or 
in  the  design  of  new  equipment.  The  two  classes  of  uses,  will  require  models  of 
somewhat  different  construction. 

In  the  area  of  operations,  equipment  models  can  be  \ised  to  optimize  the  op¬ 
eration  of  a  process.  Typically,  the  model  would  be  used  to  find  a  set  of  process 
paramters  which  result  in  the  least  variation  of  the  outputs.  Such  optimization 
procedures  might  be  used  in  lieu  of  experimental  optimization,  but  would  more 
typically  be  used  to  focus  attention  on  the  most  suitable  range  of  operation  with 
final  adjustments  performed  experimentally.  If  the  equipment  model  includes  pre¬ 
dictions  of  the  effective  disturbances,  the  optimization  procedure  using  the  model 
can  include  minimization  of  the  effect  of  the  disturbances.  This  is  an  especially 
valuable  approach  to  optimization,  since  including  disturbance  factors  in  an  exper- 
imentai  program  would  involve  an  enormous  number  of  experimental  points. 

Equipment  models  can  also  be  used  for  on-line  quality  contol  or  process  control. 
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Having  selected  the  operation  point  in  the  optimization  procedure,  the  model  can 
now  serve  to  guide  adjustments  made  locally  around  the  operating  point.  For 
example,  as  an  LPCVD  tube  drifts  with  build  up  from  successive  runs,  the  model 
cam  be  used  to  predict  the  changes  in  gas  flow  rates  needed  to  bring  the  results  back 
close  as  possible  to  the  target. 

A  model  which  captxires  the  efl'ect  of  internal  parameters  such  as  geometric 
dimensions,  and  choice  of  materials,  can  be  used  as  a  simulation  tool  for  the  design 
of  new  equipment.  Such  a  model  can  substantially  reduce  the  development  time  for 
new  equipment. 

2.3  Construction  of  Equipment  Models 

Equipment  models  may  be  constructed  by  two  distinctly  different  approaches:  phys¬ 
ically  based  mechanistic  modeling  and  statistical  modeling.  Each  approach  has  its 
distinct  advantages.  Physically  based  models  have  the  advantage  of  broad  applica^ 
bility,  good  extrapolation  beyond  the  range  of  experimental  verification,  and  good 
prediction  of  process  sensitivities.  Statistical  modeling  has  the  advantages  of  ease 
of  application  and  good  absolute  acctiracy  within  the  range  of  measmement. 

Physically  based  models  may  be  either  closed  form  or  numerical  (finite  element 
methods,  boundry  element  methods)  in  nature.  Statistical  models  are  most  effec¬ 
tively  developed  using  techniques  of  statistical  experimental  design,  such  as  Box 
“factorizd  experimental  design  and  response  surfaces"  [9]  and  Taguchi  “orthoganal 
array”  [8].  The  unifying  feature  of  designed  experiments  is  that  all  the  parameters 
of  interest  2U'e  varied  simultaneotjsly,  in  contrast  to  the  more  conventional  one  vari¬ 
able  at  at  time  experimental  techniques.  In  this  way,  the  total  experimental  range 
is  explored  with  a  minimum  number  of  experiments. 


2.4  Goal  of  the  Current  Work 

In  current  practice,  the  two  approaches  to  model  construction  discussed  above  have 
been  followed  independently.  The  broad  goal  of  this  work  is  to  fuse  the  two  methods 
and  gain  the  benefits  of  both. 

In  the  current  work,  we  have  constructed  a  physically  based  model  using  finite 
difference  niimerical  methods.  The  model  has  four  adjustable  coefficients  embedded 
in  it  which  represent  areas  of  uncertzdnty  about  the  physics.  These  coefficients  have 
been  calibrated  using  a  series  of  statistically  designed  experiments  in  order  tc  assure 
wide  range  of  application  of  the  model  with  the  minimum  number  of  expeiimenta. 
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3  MODEL  DEVELOPMENT 

3.1  Modeling  Approach 

The  model  consists  of  a  one-dimensional  representation  of  the  LPCVD  reactor  for 
polysilicon.  The  inputs  to  the  model  include:  the  physical  dimensions  of  the  pro¬ 
cess  area  (i.e.  liner  diameter,  process  tube  length,  injector  diameter);  diameter 
of  the  process  wafers;  position  of  the  W2ifer  load;  flow  rates  from  the  three  injec¬ 
tors;  positions  of  two  of  the  injectors;  temperatures  and  positions  for  each  of  the 
profile  thermocouples.  The  axial  temperature  profile  is  determined  by  a  linear  in¬ 
terpolation  of  the  temperatures  for  each  of  the  thermocouple  sites.  The  axial  flow 
ifl  modeled  with  a  convection-diff’tision  representations  incorporating  area  changes 
due  wafer  load  and  a  laminar,  plug  flow  velocity  distribution. 

The  model  predicts  the  axial  deposition  profile  of  the  polysilicon.  This  solution 
is  calculated  using  a  Newton-Raphsoa  method  on  the  center-difference  numerical 
representation  of  the  one-dimensional  system. 

The  following  sections  describe  the  analysis  behind  the  assumptions  in  the  model 
and  the  manner  in  which  it  was  constructed. 

3.2  Process  Physics 

The  process  physics  can  be  likened  to  a  coupling  between  a  number  of  physical 
mechanisms: 

•  introduction  of  gas  into  the  furnace 

•  mass  transport  in  the  axial  flow  direction 

•  mass  transport  between  the  wafers 

•  heat  transfer 


•  chemical  reaction  at  the  hot  surfaces 
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3.2.1  Injected  Gas  Dynamics 


The  injector  exit  gas  flow  is  a  highly  complex  problem.  Due  to  the  small  injector 
diameter  (-^  4  mm)  the  velocities  from  the  injectors  are  on  the  order  of  0.5  x  Mach  1. 
It  has  also  been  shown  (by  a  Poiseuille  flow  approximation  in  the  injector)  that 
there  is  a  large  pressure  drop  at  the  injector  exit.  This  indicates  that  there  is 
an  extreme  expansion  wave  at  the  exit  of  the  injector  [lO].  From  analysis  of  the 
Reynolds  number  for  the  exit  jet,  the  flow  was  shown  to  be  turbulent  (Re  ~  37). 
The  physical  understanding  of  these  systems  is  limited. 

3.2.2  Axial  Mass  Transport 

The  relative  importance  of  the  convective  flux  to  the  diffusive  flux  is  captured  by 
the  Peclet  number,  Pco  =  Where  V  is  the  estimated  average  gas  velocity  in 
the  in  the  annular  region  (V  ~  200  cm/sec),  d  is  a  relative  length  on  the  order  of 
the  furnace  length  (d  200  cm),  2md  D  is  the  diffusion  coefncient  at  the  operating 
temperature  smd  pressure  (D  ~  6000  cm^/sec).  Since  this  analysis  shows  the  Peclet 
number  on  be  of  the  order  of  one,  both  the  convective  and  diffusive  fluxes  must  be 
represented  in  the  model. 

An  indication  of  the  flow  regime  can  be  determined  by  the  the  Reynolds  num¬ 
ber,  Re  =  Where  V  is  the  estimated  average  gas  velocity  in  the  annular  re¬ 
gion  {V  ~  200  cm/sec),  is  the  gas  density  (p  2  x  10~*  kg/m*),  d  is  a  relative 
length  on  the  order  of  the  furnace  length  (d  200  cm),  and  /i  is  the  gas  viscosity 
{n  ^  2  X  10~®  kg/m  sec).  The  Reynolds  number  was  found  to  be  on  the  order  of 
one.  This  indicates  that  a  laminar  low  app»‘c‘'imation  can  be  made. 

An  analysis  was  conducted  to  estimate  the  effect  of  the  volumetric  increase  of 
gas  due  to  the  chemical  reactions  at  the  hot  surfzbces.  The  analysis  consisted  of  a 
hybrid  Poiseulle  flow  approximation  in  the  annulair  region  accounting  for  possible 
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negligible  pressure  drop  down  the  length  of  the  furnace.  The  pressure  drop  was  on 


the  order  of  0.1%  of  the  total  pressure  (~  250  mtorr).  This  means  that  a  constant 
pressure  assumption  in  the  tube  is  valid. 

3.2.3  Radial  Mass  Transport 

The  Peclet  number  can  be  used  to  indicate  the  relative  importance  of  the  convective 
and  diffusive  terms  in  the  radial  direction  by  changing  the  relative  order  of  the  length 
to  that  of  the  liner  diameter  (~  21  cm)  and  estimating  a  radial  velocity  between 
the  wafers.  Due  to  the  sjTnmetry  of  the  wafer  loading  in  the  furnace  any  convective 
flow  between  the  wafers  would  be  very  small.  From  this  analysis  it  can  be  seen 
that  the  Peclet  number  is  much  less  than  one  (Pe^  C  1),  indicating  that  convective 
transport  is  negligible  as  compared  to  the  diffusive  transport. 

Another  method  of  analysis  of  the  radial  mass  transport  effect  can  be  done 
through  the  Sherwood  number,  Sh  =  ^.  Where  fc  is  the  mass  transport  coefficient 
of  the  reaction  from  the  Arrhenius  reaction  dependence  (Ic  ~  8  x  10~*  cm/sec),  d 
is  a  relative  length  on  the  order  of  the  tube  diameter  (d  20  cm),  and  D  is  the 
diffusion  coefficient  (D  ~  6000  cm^/sec).  The  Sherwood  number  is  also  much  leas 
than  one  (Sh  C  l).  This  indicates  that  growth  rate  is  reaction  rate  limited,  or 
that  the  diffusion  time  for  the  gases  between  the  wafers  is  much  smaller  than  the 
reaction  time  at  the  wafer  surface,  indicating  that  the  radial  gas  concentration  can 
be  considered  uniform. 

3.2.4  Heat  Transfer 

An  order  of  magnitude  anadysb  for  the  relative  importance  of  radiative  to  convective 
heat  transfer  can  be  captured  for  small  temperature  variations  as:  Where 

<7  is  the  Stefan-Boltzmaun’s  constant  (5.67xl0“*  W/m*  ®K^),  T  is  the  operating 
temperature  (2’  ~  SOS^K),  d  is  a  relative  length  on  the  order  of  the  tube  diameter 
(d  —  .214  m),  and  k  is  the  thermal  conductivity  of  the  gas  (k  2x10"’  W/m  K). 
The  analysis  show  that  ~  42  indicating  that  radiation  is  the  dominant  mode 


of  heat  transfer.  This  allows  for  an  isothermal  assumption  in  the  radial  direction. 


S.2.5  Reaction  Kinetics 

The  reaction  kinetics  involved  in  the  pyrolysis  of  silicon  from  pure  silane  has  been 
studied  [2,3]  and  the  general  consensus  has  been  that  the  silane,  SiH^,  breaks  down 
in  the  gas  phase  as: 

Sin,{g)  ^  SiH,{g)  +  Hiig)  (1) 

k2 

The  SiHj  then  adsorbs  on  the  hot  surfaces,  silicon  is  deposited,  and  the  by-product, 
hydrogen,  desorbs: 

SiHiig)  ^ 

ki 

ks 

SiR2(a) 

k, 

jurf 

H2{a)  ^ 

kn 

The  reaction  equation  for  the  deposition  of  silicon  from  silane  describing  the  above 
chemical  reaction  equations  and  substantiated  by  experimental  evidence  was  for- 
mixlated  by  Roenigk  and  Jensen  [6]  as: 

_  k\Xsifj^C lot 

1  -f  KtXsifl^Ciot  +  Kf^yJCtot  (1  —  XsiHi) 

Where  is  the  molar  fraction  of  silaiie,  (1  -  Xsm*)  represents  the  molar  fraction 
of  hydrogen,  Ctot  is  the  total  concentration  of  the  gases,  K,  and  Kh  are  the  equi¬ 
librium  coefficients  (adsorbtion  to  desorption)  of  silane  and  hydrogen  and  indicate 


molSi 


m^sec 


(5) 


SiffjCa) 


Si{s)  +  .ffj(a) 


H2{g) 


(2) 


(3) 


(4) 


a. 


b. 


Figure  3:  Representation  of  the  source  and  center  injector  empirical  formulations: 
a.  Flat  model,  b.  Ramp  model. 

the  effect  of  silane  and  hydrogen  concentrations  on  the  silicon  deposition,  and  Jbi  is 
an  Arrhenius  reaction  rate  dependence  represented  as: 

=  (6) 

Where  is  the  Arrhenius  reaction  constant,  AEa.  is  the  activation  energy,  R  is  the 
universal  gas  constant,  and  T  is  the  temperature.  This  is  the  major  contributor  to 
the  effect  of  temperature  on  the  reaction. 

3.3  Model  Physics 

The  model  captures  the  dominant  charaicteristics  of  the  afore  mentioned  physical 
mechainisms  and  the  physical  representation  of  the  system. 

Gas  is  injected  through  three  injectors.  The  load  injector,  which  is  directed  in 
the  flow  direction,  is  modeled  by  a  constant  concentration  distribution  across  the 
tube  cross-sectional  area.  This  is  possible  because  there  are  no  wafers  in  the  load 
region  and  the  gas  will  diffuse  rapidly  in  the  furnace.  The  two  counter-flow  injectors, 
center  and  source,  are  modeled  by  incorporating  as  much  physical  knowledge  of  the 
flow  pattern  as  possible.  Two  empirical  flow  models  where  proposed:  a  flat  injector 
model.  Figure  3a,  and  a  ramp  model.  Figure  3b.  The  ramp  model  is  more  physically 
representative  because  it  incorporates  the  higher  concentration  of  the  gas  at  the 
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injector  exit  which  would  result  form  an  expansion  wave.  The  flat  injector  model 
gives  an  average  approximation  of  the  occxiring  exit  jet  phenomena.  Choice  of  the 
best  injector  spray  model  will  be  discussed  in  Section  5.2. 

Due  to  the  complexity  of  the  flow  physics  the  lenght  of  the  injector  spray  was 
left  as  an  adjustable  parameter; 


^  —  COnAQ,‘ny«e#«rf  (^) 

Where  Qinjteud  h*  flow  rate  of  the  gas  from  the  injector,  con  is  an  adjustable 
parameter,  A  is  the  wafer  spacing,  and  I  represents  the  effective  distance  of  the  gas 
spray.  A  linear  relationship  between  the  spray  length  and  the  flow  rate  was  chosen 
since  the  exact  physics  was  too  complex  to  determine  the  proportionality  and  a 
linear  relationship  is  more  stable  than  other  possible  relationships. 

Incorporation  of  the  injector  model  into  the  LPCVD  model  was  done  by  inte¬ 
grating  a  silane  generation  term,  J,  over  a  transversal  section  of  the  tube  to  give 
an  empirical,  one-dimens ionzd  injector  function: 

F{z)  =  [  TdLt  (8) 

J  Li 

where  Lt  is  the  perimeter  of  the  tube.  This  gives  a  one-dimensional  embodiment  of 
the  gas  injection  per  unit  length  of  tube  area. 

The  axial  mass  transport  In  the  furnace  was  modeled  by  a  convection-diffusion 
flow  representation  with  a  laminar,  plug  flew  velocity  distribution.  The  basis  for 
this  representation  arises  from  the  analysis  of  the  axial  Peclet  number  (Pe^  ^  1) 
and  the  Reynolds  number  (Re  1). 

The  radial  mass  transport  is  neglected  (Pe,  <  1  and  Sh  C  1).  This  allows  a 
one-dimensional  numerical  analysis.  Meaning  that  the  gas  concentration  at  each 
axial  position  in  the  furnace  has  a  constant  radial  gas  concentration. 

Since  the  heat  transfer  is  radiation  dominated,  as  previously  shown,  an  isother¬ 
mal  assumption  cam  be  made  in  the  radial  direction.  This  decouples  the  temper- 
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atxire  from  the  model  equationa  since  only  the  axial  temperature  will  effect  the 
diffusion  cuefficents  of  the  gases  and  the  gas  densities. 

The  reaction  model,  Equation  5,  is  included  as  a  silane  depletion  term  in  the 
silane  flux  equation.  The  model  is  combined  with  a  reaction  area  per  unit  length 
term  which  includes  the  depletion  mass  from  the  wafer  load  and  the  hot-wall  sur¬ 
faces.  This  reaction  area  per  unit  length  term  was  fo\md  by  integrating  the  areas 
of  the  wafer  surface  and  the  cylindrical  wall  surface: 

where  is  the  area  associated  with  the  surface  of  the  wafer  and  At  is  the  area 
associated  with  a  cylinder  having  a  length  A/2,  where  A  is  the  wafer  spacing,  and 
having  a  radius  equal  to  the  radius  of  the  tube.  The  evaluation  of  the  integrals  give: 

R={rl  +  r.  a)  — ^  (10) 

which  can  be  written  in  a  compact  form  as: 

R^LrZ  (11) 

where  Lr  is  the  effective  reaction  length  defined  m  ^  This  is  the  one¬ 

dimensional  depletion  term  xised  in  the  model. 

3.4  Mathematical  Formulation 

The  model  was  formulated  based  on  the  convective  and  diffusive  fluxes  of  the  silane 
and  hydrogen  gases  and  the  stoichiometry  of  the  chemical  reaction: 

V(Ws.H.)  =  0 

v(j9h.)  =  0 

Where  NsiH^  and  are  the  hydrogen  and  silane  molar  flows.  Proper  boundary 
conditions  sue  coupled  to  the  above  equations. 

Assuming  that  ?  is  the  coordinate  along  the  tube  and  is  an  inward  normal 
vector  to  the  generic  surface,  the  boundary  conditions  for  silane  are: 


20 


•  at  the  wafer  surfeure 


(13) 


•  at  the  hot  wall: 


=  -^  +  T{z) 


(14) 


where  the  term  JC  la  as  formulated  in  Equation  5  and  J(£)  takes  empiri¬ 
cally  into  account  that  silane  b  injected  into  tb.e  furnace  as  dbcussed  b  Sec¬ 
tion  3.2.1. 

•  at  the  load  zone: 

NsiB,-n  =  CoVo.h  (15) 

where  Vq  is  the  gas  velocity  of  the  gas  injected  by  the  load  injector  and  Cq  is 
the  molar  density  of  the  injected  gas; 

•  and  assuming  that  at  the  end  of  the  tube  there  b  no  chemical  reaction: 

^CsiBi 


dz 


=  0 


(16) 


Similarly,  the  boundary  conditions  for  hydrogen  are: 


•  at  the  wafer  surface  and  at  the  hot  wall: 


=  2je 


(n) 


where  the  2  Z  results  from  the  stoichiometry  of  the  chemical  reaction  in  which 
for  the  every  silicon  molecule  deposited  two  hydrogen  molecules  are  generated; 


•  at  the  load  zone: 


(18) 


•  and  assuming  that  at  the  end  of  the  tube  there  is  no  chemical  reaction: 


dz 
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=  0 


(19) 


> 


Next,  the  fluxes  in  Equation  12  can  be  combined  as  :  Nm  =  giving: 


V  = 


0 


(20) 


With  the  boundary  conditions  defined  by  the  linear  combination  of  the  afore  men¬ 
tioned  boundry  conditions. 

Assuming  now  that: 


^SUJa  ~  XsiH^  +  Cut 


Out 


XsiH^ 


(21) 


and 

=  —OutD^Xsi  +  Out 

where  Nut! Out  =  V  is  the  molar  average  velocity  of  the  gas  and  Cut  =  C'siH.  +  C'hi- 
The  model,  coupled  with  the  boundary  conditions  previously  introduced,  cam  be 
stated  as: 


(22) 


V  [-CutD'^XsiB,  +  Out'^X  5*^4)  =  0 

V^-CutD^XH.  +  CutyXB,)  =  0 


This  model  consists  of  a  system  of  elliptic  partial  differential  equations  which  can 
be  solved  numerically. 

Equation  23  can  be  reformulated  to  take  advantage  of  more  effiecient  numerical 
technioues.  Since; 

V  (Cut'V)  =  0  (24) 

and  with  the  proper  boxmdary  conditions: 


V  {-CutD^iXsiB,  +  Xb,)  +  Out^(XsiB,  +  Xb,))  =  V  (c,„.v)  (25) 

It  follows  that  the  solution  can  be  obtained  in  terms  of  one  of  the  two  concentrations 
and  the  molar  average  velocity.  The  problem  can  thus  be  formulated  as: 

V  (C.o,v)  =  0  (26) 
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V  {-CtoiDVXs,H,  +  Ctot^XsiB,)  -■=  0  (27) 

From  the  asdumptions  found  in  Section  3.2  the  equations  can  be  simplified  to  one- 
dimensionai  form  in  the  axial  direction.  Since,  in  the  one-dimensional  approxima¬ 
tion,  the  gas  flow  is  uniform  across  any  radial  area,  Equations  27  and  26  can  be 
stated  as: 


-  ^  +  A„C,„^Xsw.]  =  -  (R  -  f  (z)) 


=(R  +  RW)  (W) 

where  is  the  cross  sectioned  flow  area  defined  according  to  the  position  in  the 


tube  and 


-$£  =  v 


is  a  potential  flow  approximation  of  the  molar  average  velocity  used  for  numerical 
efficiency. 

The  boundary  conditions  associated  with  Equation  29  can  be  stated  as: 


“■  |»=o  ~  I  1  =  ^0 


tl’iL)  =  0  (32) 

The  first  condition  states  that  at  the  load  zone  the  total  gas  velocity  is  determined 
by  the  amount  of  silane  introduced  and  the  second  condition  sets  the  reference  value 
for  the  potential.  The  boundary  conditions  associated  with  Equation  28  are: 

~  ~  Qo  (33) 


dXsiH, 

dz 
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Which  again  state  that  the  silane  flow  at  the  load  sone  is  detennined  by  the  amount 
of  silane  introduced  through  the  load  injector  and  that  at  the  very  end  of  the  furnace 
the  silane  concentration  is  constant.  In  these  equations  L  is  the  length  of  the  tube 
and  Qo  is  the  quantity  of  SiH^  introduced  at  the  front  of  the  tube  for  each  time 
unit. 


3.5  Numerical  Technique 

The  discretization  of  the  equatioris  consisted  of  a  center  difference  method  over  a 
one  dimensional  grid  which  consisted  of  N+1  grid  points  having  coordinates  z,,  *  = 
1,  N+l.  Equations  28  and  29  were  decoupled  and  then  linearized  using  the  Newton- 
Raplison  method,  Appendix  A.  The  resulting  discretized  equations  were  arranged 
as  tri-diagnol  matrices  and  solved  \ising  a  standard  solving  package.  Included  in 
Appendix  B  is  the  code  of  the  polysilcon  model.  The  model  is  written  in  C  with 
two  Fortran  subroutines. 

3.6  Adjustable  Parameters 

There  are  four  adjustable  parameters  in  the  model  which  will  be  fit  to  experimental 
data  collected  from  a  designed  experiment. 

Three  of  the  constants  are  in  the  reaction  model  depicted  m  Equation  5.  They 
are  ki.K.  .  and  Kk.  These  parameters  are  physical  constants  which  arise  from  the 
chemical  reaction  for  the  deposition  of  polysilicon.  The  Arrhenius  type  dependency 
of  the  reaction  rate  is  found  in  ki,  where  ko  from  Equation  6  is  to  be  extracted 
from  experimental  data.  The  effects  of  the  silane  and  hydrogen  adsorption  and 
desorption  on  the  growth  rate  are  quantified  by  K,  and  Kk,  respectively. 

The  fourth  constant  is  iii  the  empirical  injector  function  for  the  spray  length, 
Equation  7. 
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Figure  4:  Taguchi  Lo  array. 

4  EXPERIMENTATION 

4.1  Motivation 

Careful  consideration  must  be  given  to  the  design  of  a  set  of  experiments.  The  pm  - 
pose  of  the  experiments  must  be  clearly  defined  so  that  useful  results  are  obtained. 
The  experiments  in  this  work  were  necessary  to  calibrate  and  test  the  numerical 
equipment  model.  Therefore,  the  design  of  the  experiments  had  to  investigate  a 
large  operating  space  so  that  the  model  would  be  as  accurate  as  possible.  A  second 
consideration  was  the  expense  of  each  experiment.  It  was  thus  necessary  to  design  a 
set  of  experiments  which  would  take  into  account  both  of  the  above  considerations. 

4.2  Experiment  Design 

A  set  of  experiments  was  designed  to  calibrate  and  test  the  accuracy  of  the  model. 
The  Taguchi  L9  orthogonal  array,  developed  by  Genechi  Taguchi  l8],  was  chosen 
for  this  experimental  design.  Figure  4.  This  design  allowed  the  investigation  of  a 
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large  operating  space,  with  four  equipment  parameters  at  three  levels,  in  only  nine 
experiments.  The  limitations  of  this  design  were  that  the  four  parameters  would 
have  to  be  independent  of  each  other  and  that  factor  interactions  could  not  be 
studied  &om  the  results.  These  are  not  inhibiting  limitations.  Independent  factors 
can  always  be  found  or  interacting  factors  can  be  combined  into  one  independent 
factor  [12]. 

The  choice  of  equipment  parameters  to  be  used  as  experiment  factors  was  abo 
important.  The  model  was  developed  to  predict  the  axial  deposition  profile.  To 
verify  the  model,  it  was  thus  necessary  to  choose  equipment  parameters  which  had 
a  dominant  effect  on  the  axial  deposition  profile. 

There  are  many  equipment  parameters  which  effect  the  axial  growth  profile. 
These  include: 


1.  pressxire  variations 

2.  flow  rates  from  each  of  the  three  injeciorb 

3.  temperature  variations 

4.  off  center  of  the  cantilever  or  liner  alignment,  which  would  create  a  static 
pressure  differential  across  the  wafers  promoting  convective  flow  between  the 
wafers 

5.  amount  of  prior  deposition  in  the  furnace 

6.  prior  condition  of  the  wafers 

7.  wafer  thickness  variations,  causing  temperature  non-uniformities  within  the 
wafers 


8.  thermocouple  aging 


9.  v/afer  spacing 
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10.  exterior  conditions,  including  humidity  and  temperature 

Before  the  equipment  parameters  to  be  \ised  as  experiment  factors  where  chosen  it 
was  necessairy  to  determine  the  effect  of  equipment  noises  on  the  deposition  profile. 
A  set  of  parameter  settings  was  found  which  gave  a  relatively  flat  profile.  These 
settings  were  the  base  line  equipment  settings.  The  furnace  was  run  with  these 
settings  for  a  number  of  runs  which  spanned  approximately  100,000  A  of  deposition 
on  the  furnace  walls.  After  a  clean  of  the  equipment,  the  base  line  wblb  run  again. 
This  set  of  experiments  indicated  that  the  repeatablity  of  the  equipment  was  not 
greatly  effected  by  equipment  noises.  Figures  of  the  measured  data  include  error 
bands  for  the  experimental  growth  rates  as  determined  from  these  experiments. 

From  these  experiments,  the  four  pariuneters  considered  to  be  the  most  domi' 
namt  as  to  affecting  the  axial  gjrowth  rate  and  allowed  ease  of  control  by  the  operator 
were:  (1)  pressure,  (2)  temperature,  (3)  gas  flow  rates  from  the  three  injectors,  and 
(4)  positions  of  the  two  moveable  injectors. 

This  gave  seven  parameters  which  could  be  varied.  The  L®  design  only  allowed 
four.  To  quantify  the  the  Arrhenius  and  equilibruim  coefficients  in  the  reaction 
model,  Equation  5,  the  pressure  was  a  necessary  parameter.  The  second  considera¬ 
tion  was  that  the  parameters  had  to  be  independent  of  each  other  and  the  physical 
limitations  of  the  equipment  had  to  be  kept  in  mind.  It  was  necessary  to  choose  a 
total  flow  rate  which  wouid  allow  the  desired  variation  in  the  pressure.  Based  on 
experience  of  the  BTU  Engineering  staff,  a  flow  rate  of  160  seem  was  chosen  as  the 
total  flow.  This  meant  that  the  three  injector  flow  rates  would  be  coupled.  The  load 
and  center  injector  flow  rates  were  chosen  as  parameters,  therby  fixing  the  source 
injector  flow  rate.  The  fourth  parameter  chosen  was  the  source  injector  position. 
Since  the  position  of  this  injector  effects  the  jocial  profile  uniformity  more  than  the 
center  injector  positioii.  The  factors  thus  chosen  were: 

•  pressure 
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•  load  Injector  flow  rate  {Qt^) 

•  center  injector  flow  rate  (Qe«ru«r) 

•  poeition  of  source  injector  (Xtouret) 

Factor  levels  for  the  experiment  factors  were  chosen  based  on  the  base  line 
settings.  The  high  and  low  values  for  the  levels  were  determined,  based  on  the 
prior  experience  of  the  BTU  Engineering  staff,  to  give  a  large  representation  of  the 
operating  space  for  the  furnace. 

An  independent  set  of  experiments  was  conducted  using  flow  through  the  load 
injector  only.  These  experiments  were  conducted  to  investigate  the  parameter  in 
the  reaction  model  with  the  maximum  amount  of  concentration  variation.  These 
experiments  are  discussed  in  Section  5.3. 

4.3  Procedure 

The  experiments  were  conducted  in  a  commercial  BTU  Engineering/Bruce  Sys¬ 
tems  7351C  horbontal,  hot-wall  reactor  with  three  zone  temperatxire  control.  The 
equipment  consisted  of  a  230/240  mm  x  88.5  inch  quartz  process  tube  fitted  with  a 
214/220  mm  x  70  inch  quartz  liner.  The  furnace  had  a  three  zone  heater  with  a  32 
inch  flat  zone.  The  total  ^fer  load  consisted  of  150,  6  inch  wafers,  in  six,  6  inch, 
50  wafer  quartz  boats  with  3/32  inches  center  to  center  spacing.  Only  25  wafers 
were  inserted  per  boat  giving  a  wafer  spacing  of  3/16  inches  center  to  center.  The 
first  and  last  boats  were  dummy  boats.  The  production  load  thus  consisted  of  100 
wafers  in  the  center  four  boats. 

The  recipe  used  to  run  the  fxirnace  was  the  recipe  which  BTU  recommends 
for  furnace  operation  when  running  a  flat  polysilicon  process.  The  experiments 
were  run  by  varying  parameters  in  the  recipe  according  to  the  Lg  experiment  design 
structure.  The  deposition  time  in  each  experiment  was  75  min.  Thirteen  test  wafers 
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were  inserted  at  locations  20,  26,  35,  45,  55,  65,  75,  85,  95,  105,  115,  124,  and  130 
in  the  wafer  load  of  150  wafers  (wafers  20  and  130  where  in  the  dummy  load). 

These  w^lfers  were  then  measured  with  a  nanospec  and  an  ellipsometer.  Mear 
surements  were  made  at  the  top,  bottom,  center,  left,  and  right  of  the  wafers.  The 
average  of  these  measurements  was  the  growth  attributed  tc  each  of  the  wafers. 

4*4  Experimental  Process  Optimization 

The  experiment  factors  from  the  L®  design  were  optimized  according  to  Taguchi’s 
signal-to-noise  ratio  (SN)  criteria  for  *^ominal  is  best.” 

SN  =  10  log  (  )  (35) 

\  variance/ 

The  SN  is  calculated  by  determining  the  mean  of  the  axial  deposition  profile  and 
calculating  the  variance  of  the  measured  data  from  this  mean.  The  "nominal  is 
best”  SN  is  used  because  the  optimum  deposition  profile  would  be  flat  with  the 
least  vzo-iance  about  the  mean.  This  criteria  for  the  SN  gives  a  measure  of  the 
relative  flatness  of  the  profile.  The  assumption  is  that  the  mean  can  be  scaled  to 
the  desired  value.  In  this  case,  the  mean  can  be  scaled  by  time  to  obtain  the  desired 
amount  of  total  depostion. 

Table  1  shows  the  calculated  signal- to-noioe  ratios  for  the  L®.  Figure  5  and  Fig¬ 
ure  6  indicate  the  profile  differences  between  the  best  factor  level  set  (Experiment 
2)  and  the  worst  factor  level  set  (Experiment  9). 

Choice  of  the  optimum  factor  level  settings  was  determined  by  averaging  the  SN 
for  the  experiments  in  which  the  particular  factor  level  was  used.  For  example,  the 
SN  for  experiments  1,  2,  and  3  were  averaged  to  obtain  the  relative  SN  attributed  to 
factor  1  level  1.  The  optimum  settings  were  the  factor  levels  with  the  largest  SN’s. 
Table  2  gives  the  optimum  paxaimeter  levels  determined  by  the  analysis  and  Figure  7 
shows  the  corresponding  deposition  pronie  using  these  settings.  This  optimum 
profile  is  more  axially  imiform  than  the  previous  L®  experiments  sind  the  base  line, 
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Figure  5:  Deposition  profit  with  largest  SN  for  Lj  array. 
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Figure  6:  Deposition  profile  with  Lowest  SN  for  Ls  array. 
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Figure  7:  Deposition  profile  for  optimum  factor  levels. 


thus  indicating  that  the  equipment  was  optimized. 
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5  CALIBRATION  OF  CONSTANTS 

6.1  Parameter  Fit  to  L9  Array 

The  kinetic  constants  in  the  reaction  rate  model,  Equation  5,  and  the  injector  con¬ 
stant,  Equation  7,  were  fit  using  a  design  optimization  package  called  OPTDES 
developed  at  Brigham  Young  University  [13].  This  package  \ised  a  non-linear  re¬ 
duced  gradient  method  to  optimize  a  least  square  objective  fimction.  The  objective 
function  included  all  nine  experiments  by  summing  the  squared  differences  between 
the  13  measured  data  values  and  the  predicted  values  at  the  corresponding  wafer 
Icactions  for  each  individusil  experiment  in  the  and  then  summing  the  values  for 
all  nine  experiments. 

The  cor.stants  obtained  for  the  reaction  model  were: 

ki  =  1.202x10^°  exp  (—18, 500/r)  mol/m’/s/atm 
K,  -  0.386xl0*atm-^ 

K^^  =  1.904x10^  atm" 

The  kinetic  constants  found  from  this  regression  are  substantially  different  from 
those  obtained  by  Roenigk  and  Jensen  [6]: 

ki  =  (1.6  ±  0.4)xl0®  exp  (— 18, 500/T)  mol/m*/s/atm 
K,  =-  (0.7±0.1)xl0®atm-^ 

Ks  =  (0.6  ±  0.3)xl0*atm-^/* 

Our  kinetic  constants  indicate  that  hydrogen  has  a  much  greater  effect  in  inhibiting 
the  reaction  rate;  if/,  is  three  orders  of  a  magnitude  larger.  The  Arrhenixis  constant 
in  ki  is  larger  by  an  order  of  magnitude  due  to  the  large  increas  of  if/,.  We  also 
found  that  Kg  has  relatively  little  effect  on  the  predicted  profiles  in  the  range  of 
0.01  x  10*  atm"^  to  1.0  x  10®  atm"*  since  the  denominator  is  dominated  by  the 
effect  of  hyuxogen.  The  activation  energj'  in  ki,  Equation  6,  was  not  fit  to  the  data, 
but  was  determined  by  other  investigators  [14] . 
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The  raonp  injector  gave  slightly  better  maximization  of  the  objective  function 
criteria.  The  injector  fitting  constant  of  Equation  7  was  found  to  be  5.165.  To  give 
an  idea  of  what  this  means,  for  an  injector  flow  rate  of  60  std  cm’ /min  the  modeled 
injected  gas  length  is  approximatly  3cm. 

5.2  Choice  of  Injector  Model 

The  choice  of  the  model  to  represent  the  gas  spray  dynamics  from  the  two  moveable 
injectors  was  done  based  on  two  criteria:  a  qualitative  analysis  of  the  injection 
dynamics;  and  maximization  of  the  least  square  objective  function. 

After  the  the  model  parameters  were  fit  to  the  Lq  data  it  was  foxmd  that  the 
injector  constant  for  the  flat  injector  function  was  approximatly  twice  that  of  the 
ramp  function:  9.626  (flat)  v.  5.165  (ramp).  This  is  consistent  with  the  idea  that 
the  flat  model  represented  an  average  approximation  of  the  ramp  model.  The  ramp 
function  gave  a  slightly  larger  objective  function  result,  indicating  that  the  variance 
from  the  measured  data  was  less  for  this  function. 

Since  qualitatively  the  ramp  hmction  is  more  representative  of  the  expected  gas 
dynamics  and  the  objective  function  criteria  was  maximized  for  the  ramp  function, 
the  ramp  function  was  chosen  as  injector  model. 

5.3  Calibration  Check 

To  check  the  me  of  the  Lo  experiment  set  for  calibration  of  the  kinetic  paramters 
a  ‘hnini-experiment’*  matrix  was  designed  to  capture  the  main  effects  of  the  kinetic 
parameters  in  Equation  5.  The  design  consisted  of  three  experiments  in  which  the 
pressure  and  temperature  were  varied.  In  these  experiments,  the  total  flow  was 
introduced  through  the  load  injector  only.  This  was  done  in  order  to  maximize  the 
etlect  of  the  ailane  depletion  and  hydrogen  generation  on  the  deposition  profile.  The 
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Figure  8:  Second  parameter  fit  of  reaction  model  -  Experiment  l. 

constants  found  in  the  parameter  fit  to  this  data  are: 

ki  ss  1.412x10^°  exp  (—18, 500/r)  mol/m*/s/atm 
Kg  =  0.368xl0®atm”^ 

Kh  -  1.814xl0^atm-^/’ 

These  constants  are  very  close  to  those  obtained  by  the  paramter  fit  to  the  Lg  design. 
The  difference  in  the  predicted  profiles  can  be  seen  in  Figures  8  •  10,  where  the 
predicted  profiles  with  the  respective  constant  fits  are  plotted  versus  the  measured 
data  for  each  of  the  experiments  from  three  experiment  design.  This  small  variation 
in  the  two  predicted  profiles  indicates  that  the  L9  design  was  adequate  to  fit  the 
reaction  kinetic  parameters.  From  Figures  9  and  10,  which  were  run  at  the  same 
pressure,  but  different  temperatures,  it  can  be  seen  that  the  activation  energy  used 
is  representative  of  the  system. 
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Figure  9:  Second  parameter  fit  of  reaction  modei  -  Experiment  2. 
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Figure  10:  Second  parameter  fit  of  reaction  model  -  Experiment  3 
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MATfItX  1 


6  MODEL  RESULTS 

The  predicted  deposition  profiles  closely  resemble  the  measured  profiles.  Figures  11 
-  20  give  the  predicted  deposition  profiles  for  each  of  the  Lo  experiments  and  the 
optimized  parameter  run  according  to  the  previous  Taguchi  analysis  with  the  mea¬ 
sured  data.  it  is  important  to  note  that  in  each  simulation,  the  model  was 

run  with  the  same  parameter  settings  as  the  equipment. 

The  primary  ptirpose  of  our  model  is  to  predict  the  thickness  variaton  down 
the  length  of  the  tube.  Given  a  combination  of  process  parameters  (flow  rates, 
injector  positions,  pressure,  temperature,  etc.)  a  calibration  experiment  can  be 
run  to  determine  the  amo\mt  of  time  needed  to  determine  a  target  thickness.  The 
predicted  profiles  are  therefore  normalized  to  the  means  of  the  respective  runs. 
Normalization  was  achieved  by  shifting  the  entue  predicted  curve  by  a  multiplicative 
constant  equal  to  the  ratio  of  the  measured  profile  mean  to  the  predicted  profile 
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Figure  i2:  Experiment  2. 
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Figure  13;  Experiment  3 
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Figure  20:  Optimized  factor  levels. 

mean  for  each  profile.  In  the  worst  case,  Experiment  1,  the  mean  of  the  predicted 
profile  was  10%  lower  than  the  mean  of  the  experimental  data.  The  average  error 
in  the  prediction  of  the  mean  deposition  rate  was  &%.  No  systematic  explanation 
was  found  for  the  deviation  in  the  mean  growth  rates. 

The  peaks  in  the  predicted  profiles  are  due  to  the  approximation  of  the  injector 
gas  spray.  It  should  be  noted  that  the  predicted  profiles  give  an  accurate  represen¬ 
tation  of  the  measured  profiles  except  at  positions  localized  above  the  injector  exit 
positions. 

As  a  measure  of  the  accuracy  of  the  model  to  predict  the  variations  in  the 
equipment  parameters,  the  signsJ-to-noise  ratios,  according  to  Equation  35,  of  the 
measured  profiles  were  plotted  against  the  SN  of  the  predicted  deposition  profiles. 
Figure  21.  A  perfect  correspondence  would  have  resulted  in  a  line  with  a  slope  of 
one  and  a  y-interccpt  of  zero.  The  results  (slope  =  1.05,  y-;ntercept  =  -7.1)  indi¬ 
cate  that  the  model  accurately  predicts  the  effect  of  variations  in  the  equipment 
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Comparison  Batwaon  Exparimantal  and  Pradictad  slgnal*to-nolsa  ratios 


Figure  21:  Signal-to-noise  ratio  measure.  The  numbers  indicate  the  experiment 
number  from  the  L®. 
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parameters  on  the  axial  deposition  profile.  The  slope  shows  that  the  mesured  data 
and  predicted  data  correlate  very  well  since  they  follow  the  same  trend.  The  y* 
intercept  is  a  measi'jre  of  a  constant  offset  of  the  predicted  data  from  the  measured 
data.  Meaning,  in  this  case,  that  the  model  predicts  higher  SN  for  each  of  the  exper¬ 
iments  than  those  found  &om  the  measured  data.  This  plot  shows  that  the  optimum 
operating  point,  if  determined  from  the  experimental  data,  would  be  Experiment  2, 
but  if  it  was  to  be  determined  from  the  model  predictions,  Experiment  8  would 
be  chosen.  These  values  are  very  close,  there  is  not  much  deviation  &om  either  of 
these  choices  with  respect  to  the  profile  flatness.  The  experiment  optimized  by  the 
Taguchi  SN  analysis  is  closest  to  Experiment  2. 


7  CONCLUSIONS 


A  model  has  been  developed  which  predicts  the  wafer-to-wafer  deposition  rate  of 
polysilicon  down  the  length  of  a  horizontal  tube  furnace.  Inputs  to  the  model 
include:  silane  flow  rates  from  three  injectors,  two  injector  locations,  locations  of 
and  temperatures  at  three  thermocouples,  operating  pressure,  the  number  of  wafers, 
wafer  diameter,  the  location  of  the  wafer  load,  and  other  physical  dimensions  of  the 
furnace  such  as  tube  length  and  liner  diameter.  The  model  construction  consists 
of  a  one  dimensional  finite  difference  numerical  representation  of  the  convective 
and  diffusive  fitixes  of  silane  and  hydrogen.  Silane  is  injected  and  hydrogen  is 
generated  by  reaction  at  the  wafers  and  tube  wadi.  The  silane  injection  and  mixing 
was  modeled  with  a  ramp  fimction  which  is  an  approximate  model  incorporating  a 
qualitative  understanding  of  the  injection  phenomena. 

PaLrameters  in  the  reaiction  kinetics  model  and  injector  fimction  were  fit  to  a  set 
of  nine  statistically  designed  experiments  which  varied  four  parameters,  two  injector 
flow  rates,  one  injector  position,  and  pressure,  over  tliree  levels.  These  parameters 
were  fit  to  all  nine  experiments  in  the  design  to  give  the  best  possible  fit  to  all 
the  data.  The  reaction  kinetic  parameters  were  independently  fit  to  a  second  set 
of  three  experiments  in  which  flow  was  admitted  through  only  the  load  injector. 
The  experiments  were  designed  to  explore  a  different  region  of  operation,  spanning 
a  wide  range  of  silane  and  hydrogen  concentrations.  The  close  correspondence 
between  the  reaction  parameters  fit  to  each  set  of  experiments  demonstrated  that 
the  model  extrapolates  well  to  regions  beyond  the  scope  of  the  initial  experimental 
space. 

An  extension  of  the  model,  which  included  equipment  disturbances  or  noises, 
was  used  to  optimize  the  process  by  finding  the  settings  at  which  the  process  was 
mnet  rob’jst;  that  is,  the  settings  which  gave  the  flatest  profile  ever  the  full  range  of 
the  noises.  The  settings  predicted  by  the  model  correspond  very  closely  to  an  ex¬ 
perimentally  determined  robust  operating  point  as  determined  by  a  Taguchi  signal- 
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to-noise  ratio  analysis  of  the  statistically  designed  experimental  space. 

For  the  process  engineer,  the  model  is  well  suited  for  process  optimization. 
Eqmpment  and  process  disturbances  can  also  be  included  to  determine  the  most 
robust  operating  point,  as  has  been  demonstrated.  Furthermore,  the  model  can  be 
extended  for  use  in  on-line  quality  control.  Having  optimized  the  operating  point 
by  process  optimization,  the  operator  can  use  the  model  to  correct  the  equipment 
settings,  based  on  product  measurements,  to  maintain  the  required  deposition  spec¬ 
ifications. 

The  equipment  designer  will  find  the  model  useful  for  testing  new  ideas  in  equip¬ 
ment  design.  Minor  extensions  can  be  made  to  the  model,  such  as  adding  injectors, 
to  determine  the  benefits  of  these  changes  without  the  cost  of  materials  and  exper¬ 
imentation. 
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A  APPENDIX  A  -  Numerical  Technique 


In  this  Appendix,  the  numerical  techniques  for  the  solution  of  Equations  29  and 
28  are  presented  and  discussed.  It  b  assumed  that  the  equations  are  dbcretized 
over  a  one-dimensional  grid,  which  consists  of  N  +  1  grid  points  having  coodinates 
z,-,*  =  l,iV  +  1.  The  intervab  defined  by  two  grid  nodes  t  and  »  +  1  are  referred  to 
by  using  the  standard  notation. 

A.l  The  total  velocity  equation 

Equation  29  is  a  non-linear  Pobson’s  equation  which  can  be  solved  using  first-order 
finite-difference  techniques.  In  order  to  explicit  the  non-linearity  of  thb  equation 
we  can  reformulate  Equation  28  and  29  as  follows: 

•  as  a  first  step  it  b  possible  to  observe  that  if  Nstg^  CCtoiD^^,^CtotXsiE, 
then  we  have  that  XsiH^  ^  exp 

•  Generalizing  the  previous  observation,  m  general  we  can  change  the  variable 

XsiBt  variable  4>  defined  as: 

0  =  V'  +  -D  log  XsiB^  (36) 

which  gives 

XsiB^  =  (37) 

•  Using  thb  formulation,  the  recombination  term  can  be  expressed  as  a  function 
of  the  potential  V'  and  the  non-linear  problem  solved  by  using  the  Newton- 
Raphson  method. 

To  do  thb  we  have  to  compute  the  derivative  of  the  recombination  term  with 
respect  to  the  potential  0: 

dk  _  ^  CtotXsiB^  (fea  —  fea  (Giot(l  —  XsigJ)  /2^  ^ 

1  +  kiCtotXsiB.  +  k^y/Ctotil  -XsisiT)  J  ^ 
The  Newton- Raphson  method  can  be  expressed  as  a  succession  of  linear  prob¬ 
lems.  Assuming  that  rp  =  +  Axp  we  get: 

d  f  dAtp\  dR  d  ( 

-Tz  1^“  ^^-Tz  1^“ 

so 
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•  The  discretization  of  the  previous  equation  gives  for  i  =  2,  N  —  1: 


- 

{^)  U-  A*  =  (^)  + 

(».-^~)  i*”  +^> 


(40) 


•  The  boundary  conditions  can  be  taken  into  account  as: 


—  Boundary  condition  at  the  load  zone  for  t  =  1: 

(-C-)  -  (^)  U-  A,(..  = 

I*"  (z;^) 


(41) 


—  Boundary  condition  at  the  end  of  the  tube.  The  determination  of  the 
potential  at  node  i  =  N  has  to  be  done  xising  the  following  equation: 

(;^)  {^)  - 

(f )  U-  =  -AC.  (s;Tfe)  +  («) 

+  (■'*  1*° 

A. 2  The  convection-diffusion  equation 


The  discretization  of  the  total  velocity  equation  can  be  derived  in  two  different  ways: 
either  by  assuming  that  the  potential  is  linear  over  each  element  or  assuming  that 
the  velocity  is  constant  over  each  element.  While  the  two  approaches  are  equivalent 
in  the  total  velocity  equation,  the  discretization  of  the  convection-diffusion  equation 
\b  different  in  these  two  cases.  More  technically,  the  assumption  about  the  field  can 
be  stated  in  a  framework  related  to  the  mixed  finite-element  techniques.  This 
approach  is  more  xiseful  for  the  convection-diffusion  equation  since  it  allows  the 
imposition  of  the  conservation  of  the  masses  in  a  very  straightforward  way. 

Assuming  that  Nsitr.  is  constant  over  each  element;  we  have: 

+  CtotXsiH,^  =  N 
dz  dz 
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Multiplying  each  side  of  il.c  previous  equation  by  exp{ip/D)  we  obtain: 

dXsiH,  ,  ^  = 


_  d(^XsiRtd) 
CtoiD-^—r - ^  (44) 


dz 


Integrating  both  sides  of  the  previous  equation  and  taking  Into  account  that  xp  is 
linear  over  each  interval,  we  have: 

^  dz  = 

/r  ‘  CtoiD^^^i^dz 
Evaluating  the  previous  Integrals,  we  have: 


(45) 


^'=  II"‘ 


(46) 


V'i+i " 

Rewriting  the  previous  expression  and  using  for  clarity  X  for  XsiSt 
the  flux  over  each  element  can  be  expressed  as: 

^+1  ^  V  e — _  1  — -1/ 


Defining 


D  _  ~  tPi+l 

e — -  1 


we  obtain  that  the  flux  going  from  node  i  +  1  to  node  »  can  be  expressed  as: 

^tot 


iV  = 


(48) 


(49) 


2.+1  -  A 

The  determination  of  the  concentration  C,  can  be  obtained  by  imposing  that  the  flux 
entering  into  the  “cell”  surrounding  node  »  is  equal  to  the  recombination*generation 
term.  In  this  way  we  have: 


\  __  f  7A;Ci  \  Si  —  l.iX,— 1  __ 


(50) 


-(ft -/•(!,)) 

Since  the  previous  equation  is  non-linear  we  can  solve  it  iteratively  by  using  a  the 
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Newton  method.  Defining  X  =  +  AX,  we  have; 


r/JA.+.C,+,N 

B,  ,,+ 1  AX,  4- 1  -  £ ,+ 1.,  AX, 

(  iAiCj  \ 

AX,— B,*,,_iAX,-i 

[\  ) 

J  ' 

(^)  i-^o  = 

f  /m  A  r^.  \  n. 


[  /  M«+i^y+i^ 

(  JA.O.  \ 

*1+1-'. 

V*,+  i-*.-i/ 

-{Ri-Fizi)) 


The  boundary  condition  can  be  taken  into  account  aa  follows: 


(51) 


•  at  the  load  zone  for  t  =  1: 


-  + (^)  u.  AX.  = 


(— 

V  «.  + 


i±i£!i±X^ 


-{ft-F(^))  +  (^) 


(52) 


•  at  end  of  the  tube  for  i  =  N  we  have  that  the  transport  is  due  to  convective 
effects  only.  This  implies  that  the  concentration  at  nodes  N  and  TV  +  1  is  the 
same  and  more  precisely  that  the  equation  which  determines  the  concentration 
at  node  N  is: 


AX,(-y^,) 

JA.C.  > 

1  AX,_i  1 

1  \  *1  +  1— *•— I  J 

U.+i-*.'-i  J 

J 

(^)  \xo  = 

f_  l"  i^.+i  \  _  /  ^AjCj  \  I  _  .  _ 


-F{z,)) 


(53) 
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B  APPENDIX  B  -  Polysilicon  Model  Code 
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In  this  Appendix  the  code  for  the  polysilicon  model  is  presented.  The  structure 
of  the  code  is  depicted  In  Figure  B.l.  The  m&in  part  of  the  code,  including  the 
input  routines,  initisJ  guess,  empirical  setup  of  the  injectors,  and  the  setup  of  the 
tri-diagnol  matrix  is  written  in  C.  Two  subroutines  used  to  solve  a  parameter  in 
the  convection-diffusion  equations,  bemoulli.f  -  written  by  Roberto  Guerrieri,  and  a 
tri-diagnol  matrix  solver,  dgtls.f  -  obtained  from  LINPACK,  are  written  in  Fortran. 

The  code  presented  also  has  a  inner  array  adapted  to  the  main  code.  This 
allowed  for  the  optimization  of  the  model  parameters  to  all  nine  of  the  experiment 


runs. 


Figure  B.l:  Code  Structure 
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/•  Aubroutln#  for  OPTDES  •/ 

^Includ*  <tt:dio.h> 

#irclud«  <mAth.^> 

/*  dafint  th«  Qlobel  constftnta  */ 

PI  4*fttan (1 . 0) 

#d«fin4i  9q(9c)  x*x 

♦define  cu(x)  x*x*x 

4dftflntt  irac(x)  x  -  floor (x) 

#d«£in«  CC  8.20f><i-5  /•  iftol«*K  */ 


/*  This  program  la  tha  com  of  a  simulator  for  tho  dapostion  of  */ 

/*  polysilicon  in  an  annular  diffusion  :^urnaca  */ 

main(argc,  argv) 
int  argc; 
char  *arvvi’}; 

( 

/*  dafina  variablaa  to  ba  uaad  •/ 
riZX  ^8pr  ^fcloaaO; 

doubla  ro(301],  ql^  qc,  qs,  t[301]/  p,  pci^  psi,  pfiratw,  11,  lid; 
doubla  fl,  %#dia,  idia,  D(301],  kl[301],  Ka[301],  Kh[301]«  C  tot(301]; 
doubla  r(3Cl},  dalw,  dais,  wafarp,  aw,  a> ,  cat,  sat,  lj:t30lT/  Atw[301]; 
doubla  X[3011 ,  ai(301],  dal_si[301],  dal_X(301],  aax_dal_ai,  max_dal_X; 
douMa  Ipt,  tl,  opt,  tc,  apt,  ta;  "  *  ** 

douola  ).inatica{),  injactionO,  Init^guaas  { ) ,  valocityO,  cony^diffO; 
double  fila^outO,  dgtsl^i),  dberl_(T;  * 

doubla  kk0,*ka0,  khO,  Ea7  con;  ** 
doubla  tp(4],  tql(4],  tqc[4],  tp8i[41,  ancount; 
char  input [100]; 

int  N,  nw,  n,  wafer  [150],  csntcr^in jact_«l«ja,  aourc-^^ln ject^alam,  i,  outer; 
int  count,  11 S;  **  •  — 

static  int  19 [10] [5]  -  ( 


(0  ,0, 

0, 

0, 

o; 

(0,  1, 

1, 

1, 

1) 

{0,  1, 

2, 

2, 

2) 

<0,  1, 

3, 

3, 

3) 

(0,  2. 

1, 

2, 

3) 

{0,  2, 

2, 

3, 

1) 

iO,  2, 

3, 

1, 

2) 

(0,  3, 

1, 

3, 

2) 
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(0,  3,  2,  I,  3), 

(0,  3,  3,  2.  1), 

>/ 

ap  •  fop«n  ("an^out . dat**,  *w*) ; 
fcloaa (ap) ; 

/*  nvs&bar  of  alazoants  in  tuba  M  */ 

N  -  299; 

/*  input  tha  nacaasary  paramatars  */ 

/*  prlntf  ( '"Input  tha  wafar  diamatar  (am):  ”) ;  */ 
gat a (input) ; 

%fdia  •  atof  (input)  /  1000.0;  /•  matara  •/ 

/*  printf  ("Input  tha  linar  lunar  dianatai  (tan):  ") ;  */ 
gat a (input) ; 

lid  -  atof (input)  /  1000.0;  /*  aatars  ♦/ 

Z*'  printf  ("Input  tha  total  langth  of  tha  furnaca  (cm):  ");  */ 

gat a (input) ; 

£1  •  atof  (input)  /  100.0;  /*  Mtara  */ 

dalz  •  fl  /  (N  >4^  X) ;  /*  dittanca  batwaan  nodaa  */ 

/*  printf ("Input  tha  poatition  of  tha  load  tharmocoupla  from  tha  load  and  (cm):  ") 
gats (input) ; 

Ipt  •  atof (input)  /  100.0;  /•  matara  */ 

/*  printf  ("Input  tha  load  control  tacparatura  (C) :  ");  */ 

gat a (input) ; 

tl  •  atof (input)  ♦  273.15;  /♦  dagraaa  Kalvin  */ 

/"  printf  ("Input  tha  poaition  of  tha  cantar  tharmocoupla  from  tha  load  and  (cm):  * 
gats (input); 

cpt  •  atof (input)  /  100.0;  /•  matara  •/ 

I 

/*  print f ("Input  tha  cantar  control  tamparatura  (C) :  ") ;  "/  | 

gats (input);  | 

tc  •  atof (input)  ♦  273.15;  /*  dagraaa  Ralvin  "/  | 

/"  printf  ("Input  tha  position  of  tha  source  tharmocoupla  from  tha  load  and  (cm):  " 
gat a (input) ; 

apt  •  atof (input)  /  I 00.0;  /•  matara  "/ 

/•  printf ("Inpvt  tha  aourca  control  tamparatura  (C) :  ");  "/ 
gats (input) ; 

ta  •  atof  (input)  273.15;  /•  dgraaa  Kalvin  ♦/ 

/"  datarmina  tha  ttl]  from  a  *300  C/m  tao^rature  loss  alopa  */ 

if  (rint (Ipt/dalz)  —  0)  ( 

Ipt  •  Ipt  ♦  dalz; 

> 

for  (1  •  1;  i<-rint  (Ipt/dalr) ;  !•►♦)  ( 
if  (  rint (Ipt/dalz)  —  1)  ( 
t(i]  •  tl; 

) 

•laa  ( 

t(il  •  tl  ♦  ( (-300. "Ipt)  /  (rint  (Ipt/dalz)  -  1)  )♦  (rint  (Ipt/dalz)  -  i) ; 

> 

) 

vaiiiptsi.aw batvrccn  tl  end  tc  "/ 
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for  (i»rlnt (lpt/dol«) ;  l<-rint (cpt/doii) ; 1^^)  { 

t(i]  •  tl  ♦  < (tc-tl) / (rint (cpt/doli) -rint (Ipt/dol*) ) ) * (i  -  rint (Ipt/dolz) ' ; 

) 

/*  tonporatur*  batwaan  tc  and  ts  */ 

for  (i-rint  (cpt/dalr)  ;l<-rint  (apt/dalr)  ( 

t[il  -  tc  ♦  {(ta*tc)/(rint(apt/dalx)-rlnt(cpt/dalt) ))*(i  •  rint (cpt/dalx) ) ; 

) 

/*  tao^aratura  batwaan  ta  and  tha  aourca  and  at  ->80  C/m  */ 

for  (i-rint (apt/dalx) ;  i<-(N>l);  !♦♦)  ( 
if  (  Tint (apt/dalx)  —  (N+I)  )  ( 
t(i]  •  ta; 

) 

alaa  ( 

t(ii  •  ta  ♦  (-00*  (fl-apt) /( (N+1) Tint (apt/dalx) ))• (i  -  rint (apt/dalx) ) ; 

) 

) 

/•  •nd  tavparatura  profiling  ■■■■■»  ■  */ 

/*  printf ("Input  tha  total  procaaa  praaaura  (torr) :  *) ;  */ 
gats (input) ; 

p  •  atof (input) /760 ,0;  /*  ctm  */ 

/*  calculata  tha  ailana  danaity  at  praaaura  and  tamparatura  */ 
for  (l-l;i<-(N^l) ;i+^)  { 

roli]  •  p  /  (CC  •  t (i] ) ;  /*  gm  mole  /  m^3  */ 

} 


/*  printf ("Input  tha  inner  dianatar  of  tha  injector  (mn)  :  ") ;  */ 
gat  a  (input) ; 

idia  •  atof (input)  /  1000.0;  /*  netera  */ 

/"  printf ("Input  tha  cantor  injactor  poaition  from  load  and  (cm):  *) ;  "/ 
gat a (input) ; 

pci  •  atof  (input)  /  100.0;  /*  mtara  */ 

/*  printf ("Input  the  aourca  injactor  poaition  from  load  and  (cm):  ");  "/ 
get a (input) ; 

pai  *  atof (input)  /  100.0;  /"  meters  "/ 

/*  printf ("Input  the  number  of  wafara  in  tha  load:  *);  */ 
gats (input) ; 

nw  -  atoi(input);  /"  no  units  */ 

/"  printf  ("Input  the  position  of  tha  first  wafer  from  tha  load  and  (cm):  *■);  */ 
get 3 (input) ; 

pfiratw  •  atof  (input)  /  100.0;  /"  akatara  "/ 

/*  printf ("Input  tha  length  of  tha  wafer  load  including  boats  (cm):  ") ;  */ 
gets (input) ; 

11  •  atof (input)  /  100.0;  /*  meters  */ 

/*  printf ("Input  tha  load  injector  flow  rate  (accm) :  •);  ♦/ 
get a (input) ; 

ql  -  atof (input)*(1.0/p)"(t(l)/(25>273.15))"(1.0/(cu(lC0.0))) 

"(1.0/60.0)*ro(l]; 

/"  gm  mole  /  sac  */ 

/"  printf ("Input  tha  center  injector  flow  rate  (seem):  ”);  "/ 
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guts (Input) ;  . 

1  ••  rint  (pci/d«lx ) ; 

qc  •  «tof (input)*(1.0/p)*(t(i)/(25^273.15))*(1.0/(cu(100.0))) 

*(1.0/60.0)TO(1)/ 

/♦  gm  mold  /  sec  •/ 

/*  printf ("Input  the  source  injector  flo%t  rate  (seem):  ■);  ♦/ 
gets (input) ; 
i  •  Tint (psi/delx) : 

qs  •  stof (input )* (1 . 0/p) • (t (i] / (25^273 . 15) )• (X . 0/ (cu (100 . 0) ) ) 

•(1.0/€0.0)*ro(i); 

/•  gm  mole  /  sec  •/ 


/*  printf ("input  kO:  •);  */ 
gets (input) ; 

kkO  •  etof (input)  ♦  le9; 

/*  printf ("input  ks:  ");  */ 
gets (input) ; 

ksO  -  stof (input)  •  le5; 

/♦  printf ("input  kh:  ");  */ 
gets  (ii.puw;  ; 

khO  -  stof (input)  "  le2; 

/•  printf ("input  Ea:  ");  */ 
gets (input) ; 

Ea  •  stof (input); 

/*  printf ("input  data  file  f3ag:  ■);  •/ 
gets (input) ; 
outer  •  atoi(input); 

/•  printf ("input  con:  •) ;  •/ 
gets (input) ; 
con  -  stof (input); 

/*  that  is  the  end  of  the  input  */ 

/*  set  taguchi  parameters  */ 


tp[l)  -  ,200; 
tp[2]  -  .25t); 
tp(3]  -  .3S0; 

tqldl  -  30.0; 
tql(2]  -  45.0; 
tql(31  -  60.0; 

tqcll)  •  40.0; 
tqc(2)  -  55.0; 
tqcI3)  -  70.0; 

tpsi(l)  -  113.3475; 
tpsi(2]  -  110,4275; 
tpsi(3)  -  123.5075; 

/*  begin  taguchi  19  array  */ 

for  (119  •  1;  119  <-  9;  119^-*‘)  { 

/*  set  the  correct  run  parameters  •/ 

outer  •  119; 

p  •  tpll?  !  119]  ( 1]  2 .  0: 
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/•  gm  s>ol«  /  •/ 


for  (i-X;i<-(N^l)  ( 

roll]  -  p  /  (CC  •  tlD); 

) 

ql  -  tqill9lll93  12]  ]  *  (1 .0/p)  •  (t  (1]  /  (25-^273.  IS) )  •  (l  .0/ (cu (100 .0) ) ) 

*(1.0/60.0)*ro[l]  ; 

/•  qm  »olo  /  890  •/ 


1  -  rlnt (pcl/d«lx) ; 

qc  -  tqc(19(119] ( 3} ) •  (1 .0/p) • (t li) / (25^273 . 15) ) * (1 .0/ (cu (100 . 0) ) ) 

•(1.0/60.0)*ro(l); 

/*  gm  molt  /  8*c  ‘/ 


psi  tp8iU9(119]  [4]]/100.0; 
i  rint  (psi/^lr) ; 

qe  -  (150.0-tql(19lll9]  (2n*tqc(19[119]  (3]))*(1.0/p)* 

(t (1) /(25>273. 15) )*(i.0/(cu (100.0))) 

•(1.0/60.0)*ro[l]; 

/*  gm  mol#  /  8«c  */ 

/*  calcuXato  nocotsary  arofts  and  raaction  langth,  all  units  in  matars  */ 

dslw  •  11/  nw;  /♦  artificial  wafar  spacing  */ 

aw  «  PI  •  sq(wdla  /  2,0);  /•  wafar  araa  */ 

ai  •  PI  •  aqddia  /  2.0);  /*  injact  araa  */ 

cat  •  PI  *  sqdld  /  2.0);  /*  crois  area  of  tuba  •/ 

sat  •  PI  *  (lid  )  *  fl;  /*  surfaca  araa  of  tuba  */ 


for  (i-1;  i<-(N+l);  !♦>)  { 

if  (  1  <  floor  (pflrstw/dalz)  II  i  >  floor  ( (pflrstw^^ll) /dalz)  )  ( 

Lr(l]  *  (2 . 0*PI/dalw) * ( (lid/2 .0) *dalw) ;  /•  raaction  langth  outslda 

wafar  ragion  */ 

) 

alaa 


brdl 

} 


(2.0«PI/dalw)*(sq(wcUa/2.0)^(lld/2.0)*dalw);  /«  raaction  langth 

in  wafar  ragion  */ 


/*  calculata  tha  Atw  for  aach  annular  araa  */ 


for  ( 1-1  ;i<-(N^l)  ;!-»-►)  ( 

if  (  1  <  floor (pflrstw/dalz)  II  1  >  floor ( (pflrstw^ll) /dalz)  )  ( 
Atw(i)  •  cat;  /♦  out  of  waffar  ragion  •/ 

) 

alsa 

Atw[i]  •  cat  -  aw;  /*  in  wafar  ragion  •/ 

) 


/*  and  of  araa  and  raaction  langth  calculations  */ 

/*  calculata  tha  positions,  from  tha  load  and,  injactora  */ 

/•  injactor  positions  by  alomant  */ 
cantar^injact^alam  -  rint(  pci  /  dalz); 
sourca^^ln  jact^alam  -  rint  (  psi  /  dalz); 

/*  and  of  positioning  •/ 

/*  Calculata  tha  constants  for  tha  runs  */ 

Klnat  'cs(N,  p,  t,  D,  kl,  Kj,  Kh,  C  tot,  kkO,  ksO,  khO,  Ea) ; 
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/*  CaII  function  for  the  einpericAl  injection  F(£)  */ 

injection (F,  «i,  dtlw,  qc,  qs,  center_in ject_elem,  »ource_inject_elem, 
M,  delz,  Atw,  p,  t,  con) ; 

/*  Cell  function  for  the  initial  guezses  */ 

init  gue$s(ql,  qc,  qs,  Lr,  11/  N,  ai,  X,  T,  delX/  Atw,  C_tot/ 

*  Kl/  Ka,  Kh/  0) ; 

/*  Main  body  of  program  including  the  loop  for  the  determination  of  the 
values  X  and  ai  */ 

max^del^X  -  200.; 
ttax^del^si  —  200.; 

while  (max  del  ai  >  0.0001*silN)  II  max_del_X  >  0 .0001*C_tot  (I) )  { 


velocity (N/ 
D/ 

si/  ql/ 
C_tot )  ; 

Atw,  X,  F,  dels. 

Lr, 

kl. 

Ka, 

Kh 

conv^diff  (N, 

ai,  ql. 

Atw,  X,  F,  delr. 

tr. 

kl. 

Ka, 

Kh 

0/  C_tot); 


max_^dei_X  -  fabs(del  X(l]); 
inax_del"si  -  faba  (dcT^ai  ( 11)  ; 

for  (i-2;  i<»(N-H);  ( 

if  (faba  (del^ai  (il )  >  tnax_del_8i)  { 
max_del^ai*-  faba  (del^ai  (ilT; 

)  ^ 

if  (faba(del  Xti])  >  max  del  X)  { 
max  del  X  -  fabstdel  xTi]); 

) 

) 

file_out(M/  C_tot/  X/  kl/  KS/  ICh,  dels,  delw,  pfiratw,  nw,  outer); 

/*  printfC  TEATS  ALI.  FOLKS  \n\n**)  ;  */ 

/*  exit(2);  •/ 

}  /•  19  loop  */ 

sp  •  fopen  ('’sn^out .  dat"/  "r")  ; 
ancount  0.0; 

for  (119  -  1;  119<-9;  119++)  { 
f  acanf  (sp,  "llOs**/  input); 
ancount  ~  ancount  atof  (input); 

) 

printfCan  -  %f\n"/  ancount); 

)  /♦  end  brac)cet  */ 

/*  in  this  subroutine  the  values  for  the  reaction  equation  and 
diffusion  constant  */ 

double  kinetics (n,  p,  t,  V,  XI,  Kb,  Kh,  C_tot/  kkO,  kaO/  khO,  Ea) 
int  n; 

double  p/  t(301)/  D(301)/  kl(301)/  Ks(301)/  Kh(301:/  C_tot(301]; 
double  kkO/  ksO,  khO,  Ea; 


/*  p  ;  total  system  pressure.  */ 
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/•  t  :  systtm  t«fBp«r«tur#«  */ 

doubl«  tR,  kOr  V«,  Vb,  N*,  Mb,  V,  M,  con,  coni; 

U*t  i: 

/*  CalculAto  the  total  concantration  •/ 

/•  Fron  Jansen's  paper:  Kh  -  (0,6  0.3)  a  10^2  atm*-l/2  and 

Ka  •  (0.7  0.1)  X  10^5  atm**-!.  The  constant  of  the  Arenhious 

reaction  equation  is  (1.6  *♦•/-  0.4)  x  10*^9  TOl/m‘^2/t/atffl.  ■/ 

/*  Also  place  the  tems  in  correct  unite  with  concentrations  */ 

for  (i-l;i<-(n^l);i>^)  ( 

con  -  (SC  *  t  ( i  1 ; 


)c0  •  Uc0*con;  /*  units  wl^ni/qin  aiole*sec  •/ 
kill]  •  kO*  exp((-Ea)/((1.987e-3)*t(iJ)); 

As[i]  •  ks0*con;  /*  units  stole  */ 

Kh(i)  •  kh0*sqxt (con) ;  /*  units  mole  •/ 

) 


/*  these  axe  the  values  and  calculation  for  the  diffusion  equation  * 

Va  •  14.3;  /*  Atosde  voliiate  hydrogen  */ 

Vb  •  46.8;  /*  Atewaic  Volume  silane  */ 

Ha  *  2.016;  /*  Molecular  weight  of  hydrogen  */ 

Hb  •  32. 117 S;  /*  Molecular  weight  of  silane  */ 

V*  aq((pow(Va,  (double)  (1./3.))  -«■  pow(Vb,  (double)  (1./3.)))); 

M  -  sqrt(  (l.O  /  Ma)  ♦  (1.0  /  Mb)  ); 

for  (i-l;i<-(n^l)  ;i«H^)  ( 

tR  •  (((t(i)-273.15)*9.0/5.0)^32)  ♦  460.0; 

D(il  ^  .0069  •  (pow(tR,  (double)  (3./2-))  /  p)  *  <  M  /  V  ); 

/*  Convert  D  (ft*2/  hr]  to  (in*2/sec)  */ 

coni  •  (1.0/3600.0)  *  (0.0929); 

0[i]  •  D(i)  *  coni; 


/*  initialise  the  C_tot  array  */ 

for  (i»l;i<»(n-H)  ;i-M-)  ( 

C  tot(i)  -  P  /  (GC  *  t(i]); 

} 

/*  end  of  routine  */ 


double  in jection (fl,  ail,  delwl,  qcl,  qsl,  centerl,  sourcel, 
n,  dels,  Atw,  p  ,t,  con) 
int  n,  centerl,  eourcel; 

double  f 1(301],  ail,  delwl,  qcl,  qsl,  dels; 
double  Atw[301],  p,  t(301),  con; 

( 

/*  Define  local  variables  */ 
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double  Xic/  119,  coni,  fi(301],  fc(301]; 
int  bnc,  bnar  •rro?c,  mxxorm,  12r  13; 

/*  constant  function  spraad  over  distance  of  spray  */ 

/*  astiinata  injaction  length  based  on  €0  seem  covers  20  wafers  */ 

/•  at  625  C  and  .250  torr  •/ 

coni  -  (con  *  delwl)  /  (SO .  0*  <760 .0/ .250)  •  ( (625•^273 .15)  /  (25^273 . 15) )  • 

(1.0/ (cu(lOO.O)))* (1.0/60.0)* 
(<.250/760)/(GC*(625^273.15)))  /  ail); 
lie  -  coni  *  (qcl  /  ail); 
lis  -  coni  *  (qsl  /  ail); 

/*  spread  injection  over  injection  length  (lie  and  lis)  */ 

bnc  -  centerl  *  floor (lie  /  delz); 
bns  -  sourcel  floor  (lis  /  dels); 

if  (  bnc  <  1  >  ( 

errorc  -  (1  -  bnc)  ♦  1; 
bnc  -  1; 

) 

else 

errorc  •  1; 

if  (  bns  <  1  )  ( 

errors  -  (1  -  bns)  1; 
bns  -  1; 

) 

else 

errors  -  1; 

/♦  find  injection  for  center  injector  •/ 

for  (il  -  1;  il<-<n+l);  il-M^)  { 
fc[il)  -  0,0; 
fslil]  -  0.0; 

} 


12  -  1; 

for  (il  -  bnc;  il<-centerl;  il-^-**)  { 

fc(ill  -  (( (2  ./ (centerl-bnc+l> )  *qcl)  /  (centerl  •  bnc  -f  2))*  12; 
12  -  5.2  >  1; 

) 

/•  find  injection  for  source  injector  */ 

11  -  1; 

for  (12  -  bns;  i2<-sourcel;  12*^)  ( 

fs[i2]  -  (( (2 ./ (sourcel*bns*l) ) *q8l) / (sourcel  -  bns  ♦2))  *  il; 
il  -  il  +  1; 

} 

/*  find  injector  amounts  for  all  positions  of  the  furnace  */ 

for  (13  -1  ;  13<-(n^l);  i3-»”^)  { 
fl(i3)  -  fcli3]  fs(i3); 

} 

/*  relocate  the  injection  sites  from  elements  to  nodes  */ 

for  (13  -  X:  i3<-(n*l);  i3^*)  { 
fl(i3]  -  ll[i3)  /  dels; 

) 
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/*  This  function  dofinos  th«  lniti«l  gu«t»  for  tho  t9«vton-Kaphion 
•oiution  of  tho  modoX  */ 

doubln  init  guotstqX,  qc,  qs,  XI,  n,  •!,  X,  F,  <UXs,  Atw,  C  tot, 
kl,  K«,  Kh,  D) 

int  n; 

doublo  ql,  qc,  q«,  Lr[30Xl,  11,  sl[30l],  Xi3011,  r(301],  dolz,  Atw(30l]; 
doubXo  C  tot(301],  kl(301),  Kjt30X],  Kh[30X;,  0[30I1; 


/•  local  vaxiablas  •/ 
int  i; 

doubXo  X^guaaa,  V(301],  raact^rata,  A [301],  Q(301]; 

/*  Firat  calculato  tho  constant  X^allana  */ 

/•  fix  X[1  */ 

for  (i-X;  !+-*-)  { 

X[il  -  (1./2.); 

) 

X^guasa  -  (1./2,) ; 

/*  calculata  tha  valocity  •/ 

/•  Tha  constant  raaction  rata  */ 
fcr  (l-l;i<-(n+l)  ( 

raact^rata  -  (kl[i3  *X_guaaa*C_tot  [i] )  /  (l^-Ka  (i]  •X_guaBS*C_^tot  [1] 

♦Khlil*aqrt (C_tot[i]^(l.  ••  X  guaaa))) 
R[i)  •  2,r  [1]  *saact_rata;  ** 

) 

/*  baaad  on  conaarvation  of  flow  ainca  tha  danaity  is  constant  */ 

Q(l]  -  ql  ♦  r[l]*dalz; 

Vll]  -  Qll)  /  (Atw(l]*C_tottl3); 

for  (i-2;  l<-n;  i-M-)  { 

Q(il  -  Qli-l]  +  <Rti)  +  Ftl])*dalx; 

V(i]  -  Q[i]  /  (Atw[i]*C_tot(i]); 

) 

/*  calculata  ai[j  */ 

altn-fl]  -  0.0; 

Cor  (i»n;  1>-1;  i — )  ( 

allij  -  alli-H]  ♦  Vli]*dalz; 

) 


/*  This  program  avaluatea  tha  valocity  aquation  for 

tha  daposition  of  pclyailicon  in  an  annular  diffusion  furnaca*/ 

/♦cccccccecccccccccccccccccccccccccccccccccccccccccccccc*/ 

/*  Starting  tha  function  valocity ()  */ 

doubla  valocity (n,  si,  ql,  Atw,  X,  F,  dalr,  Lr,  dalasi,  kl,  Ra,  Kh, 


Dt  C_tot) 

Int  n;  "" 

doubl«  dtlx,  »1[30X1,  ql,  Atw[3Cl],  X(301],  r(301J,  Lr(30ll; 

doubl*  dl*l_*H30X),  kl(301),  IU(3011,  Kh[301],  D[30X),  C_tox(30X]; 

/♦  n:  4ot  nodes; 

•!(]:  'old'  value  of  the  potential#  which  is  actually  tha 
'guaased'  al(]  for  the  firat  iteration.  On  return  of 
'velocity O'  ail]  carries  the  new  values  of  potential. 

X[]:  is  the  previously  established  concentration  profile# 
in  tha  first  iteration  it  will  be  the  initial  quass 
values . 

Fi):  Injection  array  where  F(i]  is  the  injected  flow 
into  the  ith  element.  Note  ith  elenent  is  between 
nodes  i  and  (I'M) .  */ 


{ 

/*  Declaration  of  the  local  variables  */ 
double  Si  oldiaoi],  si_new[301)#  x_old(301]; 

double  R(?01] #dA_by_dsi(301] # react^rate#  force [301]#  conl#con2#con3#con4; 
double  lower [30x7#  upper [301] #  dia9[3011; 
double  dgtsl^O  ; 
int  i,  info#*"npl; 

/^cccccccccccccccccccccccccccccccccccccccccccccccccc*/ 

/«  Assign  si[]  to  si_old[]  for  later  use  */ 

fora-l;i<-(n^l)  { 

si  oldti]  -  si[il; 

X  old[il  -  X[i] ; 


) 


/•cccccccccccccccccccccccccccccccccccccccccccccccccc*/ 
/#  dR__^by^^dai  •  dR/dii  at  ai^old  */ 

/•  react'*rate  is  tha  reaction  rate  */ 


forU-l;  i<-(nt‘l);  ( 

react  rata  -  <kl [i] *C_tot [i] *X  old(i3)/(l.  ♦  Ks[iJ«C  tot[i]»X  old[i)  ♦ 

Kh li]  ‘sgrt  (C“tot  (i]  *  (I, -X_old  [i] ) ) ) 
R[i]  •  Lr [i] *raact_rate;  *  * 


coni  -  (C_tot (i]*X_oldIi] )*(Ka[i]-(Kh[i]*pow(C  tot (i) * (1 .-X^old [i] ) # 

(double)  -0. 5)7/2.  ); 

con2  -  (1.0  ^  Ks(iJ*C_totIi]*X(i]+Kh(i]*sqrt(C_totli]*(l.  -  X_old(i]))); 
dR_by_dsili]  -  -1 . *Lr (!] • (raact.rata  /  D[i])*(1.0  •  conl/con2) ; 


) 


/•cccccccccccccccccccccccccccccccccccccccccccccccccccc*/ 
/*  now  let's  solve  the  potential  equation  */ 

Defining  the  coefficient  matrix  */ 

/*  for  i  -  1  •/ 


lower [1]  *  0.0; 

uppertl]  -  -1.0* (2,0*Atwl2]*C  tot (2 ] /del* ) • (1 . 0/del i) ; 

cli«9tl]  -  ((2.0«Atw[2]«C_totT21/d*l*)«(1.0/deli)  -  (dR_by_d3l [ 1] ) ) ; 

/"  for  i  •  2  to  n-1  •/ 


iow«v[i]  -  -1.0»(Atwli]*C_tot(il/d*i*)‘'.1.0/d«li); 
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upp«r(i)  -  •I.O*(Atw(i-i.l)»C  totti^l]/d«lf)*(l,0/d*l*); 

diagii)  •  ( (Atw(i]  tot  (1)7  (•q(d#l«) ) )  *  (Atw(l’^X)  •C^tot  (l^l)  /  (•q(dal*) ) ) 
dR  by  dJlTi] ) ;  ^ 

) 

/•  for  1  -  n  */ 

low«r(n)  -  -X .0* (Atw(n) *C_tot (n) /dol*) ♦ <1 .0/d«l*) ; 
uppor  [  n )  «•  0 . 0 ; 

dl4g(n)  •  ((Atw(n)*C  totIn]/(iq(d#Xt)))  ♦  (Atw(n-H) *C_tot (n+X) / (gqcdolx) )  ) 
dR_by_d*iTn) ) ; 

/*  Dofining  tho  forcing  function  voctor  */ 

/•  for  1  -  1  ♦/ 

conX  -  2.0»At»#(2)*C_totl2)/<loXx; 
con2  •0.0;  ^ 

con3  •  RIX)  ♦  r(l); 
con4  •  X.O  /  doit; 

forco(l)  •  conl»con4*  (si_old(2)  •  tl^oldd))  ♦  con3  +  2*con4*ql; 

/•  for  1  •  2  to  n-1  •/ 


for (i-2;i<»(n-l) ;!♦♦)  ( 

coni  *1.0  /  dolt; 

con2  •  (Atw(l)»C  tot (1) /dolt) ; 

COn3  •  R(i)  ♦  r(I); 

con4  -  (Atw(i+I)*C_tot(i^X)/dol*) ; 

forco(i)  -  -1.0«con2^conl* (•!  oXd(i)  -  mi  old[i-X])  ♦ 

{1.0*conX«con4)M*i_cld(i+X]‘*-  •!  old(i))^  con3; 

) 


/•  for  1  -  n  •/ 


coni  •  <1.0/d«l*); 

ccn2  •  (Atw(n) •C^tot (n) /doXt) ; 

con3  •  R{n)  ♦  F (n) ; 

con4  •  <Atw(n^X)*C  tot (n^l) /dolt ) ; 

forco[n)  -  -l.^conT*conl*(ai_old(n)  -  rl^oldtn*!))  •  con4*conX* (ri^oldln} ) 
♦  con3;  ^  ^ 


/*  How  tho  LIHPACK  subroutlno  la  callod  «/ 

for  (i-X;  i<*n;l*M')  ( 
lowor(i*l)  •  loworli); 
dl4g(i-X)  •  dlag(i); 
upporIi“l)  -upp«r(i); 
forcali^l)  -  forca(l]; 

} 


low€r(n)  •  0.0; 
dlagln]  •  0.0; 

uppar(n]  •  0.0; 
forcain)  -  0.0; 

npl  “  n; 

dgtal_(4nplr  lower/  dlag^  upper#  force#  ilnfo); 

/*  check  if  aolutlon  found  */ 

If  (  info  f-  C  )  { 

printf(**\n  Matrix  la  alnguler#  Info  -  %d\n"#info); 
exit (3) ; 
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/*  Atiign  th«  r«turr.«d  vaIuas  of  OGTSL  to  d«l_ai 
and  find  si\i)  at  tha  naw  time  step  */ 

for  { 

del  si(i)  •  forceli-l]; 

il_new[i]  •  si_old[i)  del_si[i]; 

•iji]  •  sl_new[i];  * 

/*  printf Csil^d)  %e  X(%dl  leXn-, i, si (i] , i, X(i] ) ;  •/ 
) 

sKn+l]  •  0.0; 


/•cccccccccccccccccccccccccccccccccccccccccccccccccccc*/ 

) 


/*  This  program  evaluates  the  Convection-Diffusion  equation  for 
the  depojition  of  polysilicon  in  an  annular  diffusion  furnace*/ 

/•cccccccccccccccccccccccccccccccccccccccccccccccccccccc*/ 

/•  Latest  version  as  of  DEC. 17th  1987.  */ 


/*  Starting  the  function  conv_di£fO  •/ 

/♦  We  need  to  specify  the  values  'si'  for  this  function  */ 


double  conv^diff(n,  si, 
Xn ,  D  / 


int  n; 


ql,  Atw, 
C  tot) 


X,  T,  dels,  Lr,  del_X,  kl,  Ks, 


double  •slc,  sl(3013,  ql,  Atv[3Clj, 
double  del_X(301J,  kl[301),  Ks[301] 


iwWAj,  t 

,  Kht301],  D(3013,  C_totl301] 


/*  n:  #of  ele  ^ats;  #  of  nodes  is  (n*^l)  going  from  0  to  n; 

si[]:  potential  calculated  in  the  function  ' velocity' for 
the  current  time  step; 

Xt]:  'old'  concentration  of  Silane,  which  is  actually  the 
'guessed'  X(3  for  the  first  iteration.  On  return  of 
' conv^dif £ 0 ' ,  Xt]  carries  the  'new'  values  of  Silane 
concentration . 

Ft]:  Injection  array  where  Fti]  is  the  Injected  flow 

into  the  ith  element .  Hote  ith  element  is  between 
nodes  i  and  (i^l) .  */ 


( 

/*  Declaration  of  the  local  variables  */ 

double  X  old(301]  X_new(3Cl],  X_tot[301],  si  old(301]; 

double  bI  lpl(30i:,  Bipl_i(301],  Bi_iml  [301]  ,*Bijnl_l  (301]  ; 

double  R(To1] ,dR_by_dX(301] ,react_rate,  £orce(30i;7 

double  con. coni, con2, con3, con4; 

double  lower [301],  upper [301],  diag[301]; 

double  dgtsl^O,  dberl_(),  bp,  bn; 

int  i,  info;* 

/♦cccccccccccccccccccccccccccccccccccccccccccccccccc*/ 

/♦  Assign  X(]  to  X^oldl]  for  later  use  */ 

for  (1*1;1<»  (n***l) ;  i-****^)  { 

X  oldfi]  -  Xli]; 
sl  old[J ]  •  Si[i] ; 

} 
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/•ccccccccccccccccccccccccccecccccccccccccccccccccec*/ 
/•  dB_by_dX  -  dR/dX  •/ 

/*  xvact  race  is  ths  r«'  ction  rats  */ 


<ot(i-l;  i-M-)  ( 

raact  rata  -  (kltU'C  x.ot[l]*x  old(l))/(l.  ♦  Ka(l)*C  tot  til  *x  old(iJ  ♦ 

Kh(l)»aqrt(C~totliJ*(T.-X_old(i]) ) > ; 


Rli]  •  laC  [1  l«r««ct_r«t«; 


coni  •  {K*lil*C_tot(il)-(Kh(i]*0.5*pow(C_tottil»(l.-X_ol<l[ll)  ,  (double)  -0.5)); 


con2  -  1.  ♦  K*(i)*C_tot(i)*X_oldCi)  ♦  Kh (i)  •sqrt  (C_tot  (i)  •  (1 . -X_old(il ) ) ; 
dK_by_dX(i]  •  Xrf  (i]  *re»ct^ra’-o*  ( (1  ./X^oldli] )  •  (conl/con2) )  ; 


) 


/♦cccccccccccccccccccccccccccccccccccccccccccccccccccc*/ 
/*  Lot  US  evaluate  the  vectors  B[]  where 
Bi^lpKl]  -  B(i,l4l],  etc. 

/*  i  node  i  •  1/  Oo  Is  Introduced  */ 

con  -  <si[l)  •  si(21)/  0(1]; 
dberl  ((con,  £bp,  4bn) ; 

Bi  ipT(l)  -bp  *  0(1] ; 

Bipl_i(ll  •  bn  *  0(1] ; 

/*  0  node  1  •  siln+l]  -  0.0  */ 


con  -  (si(n]  -  si(n-l])/  Din]; 
dberl  (icon,  4bp,  4bn) ; 

Bi_UJa(n)  -  bp  *  D(n]  ; 

Bi5Il^l(nj  «  bn  *  D(nj; 

/*  For  i  •  2,  (n-1),  the  B(]'s  are  ♦/ 

tor(t-2;  i<»(n-l);  !♦«►)( 

con  -  (si(i]  •  si(l+l])/  0(i]; 
dberl_(4conf  4bp,  £bn); 

Bl  ipl(l]  -  bp  *  0(i] ; 

Blpl^ili]  -  bn  *  D(i] ; 

con  -  (8i(i]  -  si(l-l]>/  Oti]; 
dberl_(4con,  4bp,  4bn) ; 

Bi  ialdl  -  bp  *  0(i]; 

Bimini (i]  -  bn  *  D(i] ; 

) 


/•ccccccccccccccccccccccccccccccccccccccccccccnccccrccc*/ 

/*  now  let's  solve  the  rorvection-dif fusion  equation  */ 

/*  Oefininq  the  coeff ici  .■  .  matrix  «/ 

lower (1]  •  0.0; 

upperd]  -  -(2.*Atw(2]*C  tot  (2) /delz)  ♦  ( (Bi  Ipl  (1] ) /delz)  ; 
diaq(l]  -  (2.*Atw(2]*C_tot(2]/delT)*((BipO(l])/delx)+ciR_by_dXll); 

for  (i»2;i<-(n-l)  ( 

lowerdJ  -  -(Atwd]*C_tot  d]/delz)*(Dl  iml  (i] ) /delz; 
upperd]  •  -(Atw(i>l]*C  tot  d+l] /delz)  *  (Bi  ipl  (i] ) /delz; 

diagd]  -  ((Atw(i*l]*C_totd>l]/delz)*(BipI_i(i]/delz)  +  (Atw(i]*C  tot  ( i] /delz)  * 
(BiJia  id]/delz)  ‘J-  dK  by_dX(i]);  ~ 


) 


low«r  In]  -  -  (Atw  [n]  *C_^tot  (n]  /d«lt )  •  (n) )  /d«lt; 

upp«r(n]  •0.0;  *  ** 

di*g  (n]  -  ( (Atw  [n^l]  *C_tot  (n^l)  /d*ljc)  ♦  (8i_olcl(n]  /del*)  ♦ 

(Atw(n]  •C^tot  (n)  /del*)  *  /d^lz)  ♦  dR_by^dX(n] )  ; 

/*  Defining  the  forcing  function  vector  */ 

forced]  -  <2.*Atw(2}*C_tot(2)/delx)»((Bi  ipld)*x  oldI2]-Bipl  1  (1]  •x^old  (1] )  / 

d«ll)  -  (R(iT  -  F(l])  ♦  T2.*qi/<1«11)  ; 


for  (i“2 ;  i<“  (n-1 )  ;  !+■«■)  ( 

forc«til  -  ((Atw(i+l)*C_totti+l]/di«l*) 

* ( (Bi_lpl (i] [11 ‘X  oid[i) > /d«l*) 

-  (Atw(i]«C  tot(i]/d«lz) 

•<  (BiJia_i(iT»X_oldti]-Bi  lnl(i]»X_cld(l-l])/d«li)  -  (R(ll -Fli] )  )  ; 


1 

force(n]  -  -(Atw  [n-t-l]  •C_tot  [n+1) /d«li)  •  (X  old!n]*si  oldIn]/d*li)  - 

(Atw(n) *C_tot (nl /d«ii) • ( (Bl»l_i(n] •X~old(nl-ai~imi [n] *X  old ( n-1] ) /d«li) 
-  (R(n]  -  F(n]);  -  ~  _ 

for  (i-1;  i<-n;  I-*-*-)  { 
lower (i-1)  •  lower [i]; 
dieg[i-l]  •  diagd]; 
upper {i“l)  •  upper [i]; 
force (i-lj  -  forceji); 

} 


lower [n]  •  0.0; 

<lieg(n]  •  0.0; 

upper [n]  •  0.0; 
force  in]  —  0.0; 

/*  How  the  LINPACK  subroutine  is  celled  */ 
dgtsl_(£n,  lower,  dieg,  upper,  force,  iinfo) ; 
if  (  info  !•  0  )  { 

printf("\n  Matrix  is  singular,  info  -  %d\n",info); 
exit (4) ; 

} 

/*  ASbign  the  returned  values  of  DGVSL  to  del_C 
and  find  C(i)  at  the  new  time  step  */ 

for  (i-l;i<-n;i+-^)  ( 
del_X[i]  -  £orced-l]; 

X  new(i)  -  X  oldli]  +  del  X{i); 

xTi]  -  X_newTij;  * 

/*  printf ("8i[%d]  %e  X(%d]  %e\n",i,  si[i],i,  X(i]);  */ 

} 

X{n-H]  -  Xin]; 


/•cccccccccccccccccccccccccccccccccccccccccccccccccccc*/ 

) 
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/*  This  lubroutint  finds  th«  correct  %#Af«E  ^owth  rtttat  *nd  pumps  thsm 
to  tho  fils  growth. djit  */  ^ 

doubls  filo_out (n,  C^tot,  X,  kl,  Kj,  Kh,  dsls,  d«lw,  pflritw,  nw,  outor) 
.^nt  n,  nw,  outsr;  ' 

C  totI301),  X(301),  kl(301],  JU(301),  Kh[301],  dslw,  dslr  • 
douoXo  pfirstw; 


/*  Iccsl  variAblss  */ 

FILE  ♦sp,  ♦»pl; 

Int  i,  j; 

•p  •  fopsn(*opt_growth.dAt",  "w")  ; 
spl  -  fopsn  (**opt^con.dit",  "w"); 

fprintf (tpI,*.col  nods  X\n"); 

for  (i*l;  l<»n‘*-l ;  { 

fprintf (spl, »%d  \t  %o\n",i,X[i] )  ; 


wsfsr  •  pfirstw; 

3  -  1; 

whils  (j  <•  nw)  { 

1  »  floor (wafsr/dolz) ; 
if  (  i  —  0  )  ( 

1-1; 

) 

X^rim#  •  ((X(i4.li  -  X[i] ) /1 . 0)  *  <fr4c  (wsfsr/do:!  2) )  ♦  X(i); 

C_tot_pEiii»  -  -  C_tot(i))/d«,\i)Mfr»c(w«f«/<J«lx))  ♦  C_tot(l] 

"  ((klil*!]  -  kl  [1] ) /1 . 0)  *  (f  rac  (waf«E/d*lz ) ;  •*  kill]; 

Ka_priaa  “  <(Ka[l-*-l]  -  Ka  [1] ) /l.O)  *  (frac(wa£ar/rt*li) )  ♦  Ka(i]; 

Kh_prime  -  (<Kh(l»l)  -  Kh [i] ) /I .0) • (frac (wafar/dalt ) )  .  Khll); 

raact_rato  -  (kl_prljne»C_tot_Ftiin«»X_prljna)  / 

(1.  ♦  K8_prinia*C_tot _ptlma*x_priaie  +  Kh _pri»a» 
aqrt  (C_tot_pri»i»»<l.-x_prijna) ) ) .- 

r«act_rate  -  react_rata*28. 086* (l,/2.328) • (1 ./(Cu(lOO) ) ) *60.*lal0; 

/*  A/mln  •/ 

li  (  J  —  20  )  { 

fprintf (ap,  "kfXn", raact ^rata) ; 


If  (  j  —  26  )  ( 

fprintf  (ap,  "kf \n", raact_tata) ; 


if  <  j  —  35  )  { 

fprintf (ap,  "%f\n",taact  rata); 


If  (  J  —  45  )  ( 


> 
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£pr intf (sp,  ■%£\n", ; 

) 

If  (  j  —  55  )  ( 

£ptint£(*p,  ■%f\n",r««ct_>:Bt*)  ; 


if  (  1  —  65  )  ( 

fptintf (tp,  "%£\n-, t««ct_r»t*) ; 


if  (  j  —  75  )  ( 

fptintf (sp,  "»£\n», tBsct_rat«) ; 


if  (  j  ~  85  )  ( 

fptintf (sp,  ■%£\n’,t*sct_rst«) ; 


if  (  1  —  95  )  ( 

fptintf (sp,  "Sf \n*, rssct^rsts) ; 


if  (  j  —  105  )  ( 

fptintf (sp,  •%f\n-,rssct_tBtB) ; 


if  (  1  —  115  )  ( 

fptintf (sp,  ■%f\n",rssct_rst«) ; 


if  (  j  —  12«  < 

fptintf (sp,  •%£\n»,tssct_tfcte) ; 


if  (  j  —  130  )  ( 

fptintf (sp,  "%f\n«,t*sct_rsts) ; 


j  -  j  *  1; 

wsfax  •  wafer  ♦  delw; 


) 

f close (sp) ; 
fclose(spl) ; 

la»:8t  (outer)  ; 

/*  and  */ 

) 

/•  vetiance  solver  */ 

double  least (outer) 
ii\t  outer; 


( 


FILE  ”spl,  ‘spl,  *8p3,  ‘feloseO; 
double  c,  aa,  ap,  v,  id[151,  pllS),  sn,  vp, 
int  1; 

char  l.iputllOO)  ; 
if  (  outer  «■  1  )  ( 


VB 
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«pl  •  fop«n("l.diAt- ,  T*); 

if  (  out«c  ~  2  )  ( 

•pi  •  fop#n("2,dmf,  -r*); 

if  (  out«r  “  3  )  { 

•pi  -  £op«n(*3.<Uf,  "r**); 

if  (  OQt«r  —  4  >  { 

•pi  -  £op«n(" 4.dat*,  "r")  ; 

if  (  OUt«E  •-  5  )  { 

§pl  •  f op#n (*S . dAt" #  "r")  # 


if  (  cuter  —  6  )  { 

•pi  •  f open ("6 -det"# 

if  (  outer  —  7  )  ( 

•pi  •  f open ("7.dat", 

if  (  outer  —  8  )  < 

•pi  ■  fopen (*8 . det" # 

if  (  outer  “  9  )  { 

•pi  -  fopenC9.det"# 

} 


■r*); 


•r"); 


■r")  ; 


-r")  ; 


•p2  -  fopen("opt^9rowth,det*',  "r"); 
•p3  -  fopen(*^n_put.det","e’')  ; 

/*  input  the  ■•••ured  value*  */ 


for  (i-0;i<-l2;i+'^)  { 

f*can£(*plr"'10*"r input) ; 
mil]  -  etof (input); 

) 


/*  input  the  predicted  value*  */ 

for  (i-0;i<-12;i-M.)  { 

f  •canf  (•p2,  •%10*-,  input)  ; 
p(i)  -  atof (input) ; 

/*  aolve  for  average  •/ 


am  **  0  >  0 ; 

•p  •  0.0; 

/♦  input  the  meaeured  value#  */ 

for  (i-l:i<-ll;i^^)  ( 
an  *  am  mil]; 
ap  -  ap  ♦  p(i]; 

) 

am  *  am  /  11.0; 
ap  •  ap  /  11.0; 

/*  •Give  for  the  variance  •/ 


V  ••  0.0; 

for  (i-l;i<-ll;i-^-*')  ( 

V  -  V  ♦  pow ( (p(il “ffi(il) # 

> 

V  •  V  /  XI. 0; 


(double)  2.0); 


/♦  eolve  for  SH  */ 

•n  •  -10  *  loglO (v) ; 

••si  \  ; 

f close (•p2) ; 
fclo*e(#pl) ; 
f close (epB) ; 


) 
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B.2  bernoulli.f 


SDBROCTlNi:  dber  (X,BP,DBP,BN,DBN) 


THIS  ROUTINE  EVALUATES  THE  BERNOULLI  FUNCTION 
ANP  ITS  DERIVATIVE  FOR  THE  ARGUMENTS  X  AND  -X 


INPUT  VARIABLE 

X  -  ARGUMENT  OF  THE  BERNOULLI  FUNCTION 

OUTPUT  VARIABLES: 

BP  -  BERNOULLI  FUNCTION  OF  X 

BN  -  BERNOULLI  FUNCTION  OF  -X 

DBP  -  DERIVATIVE  OF  BERNOULLI  FUNCTION  OF  X 

DBN  -  DERIVATIVE  OF  BERNOULLI  FUNCTION  OF  -X 

DATA: 

XLIM  •  iO.»*(-ND) 

EPS  -  10.**(-NDIGIT) 

ND  -  NUMBER  OF  ADMITTED  ERRONEOUS  DIGITS 

NDICIT  -  NUMBER  OF  SIGNIFICANT  DIGITS 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 
DATA  XL1M,EPS  /O.Ol,  i.D-15/ 


Executable  code. 


AX-DABS (X) 

COMPUTE  THE  ASYMPTOTIC  VALUES  OF  THE  BERNOULLI  FUNCTION 

IF(AX.GT.80.)  THEN 
IF(X.GT.O.)  THEN 
BP-0. 

DBP-0. 

BN-X 
DBN— 1. 

RETURN 

ELSE 

BP— X 
DBP— 1. 

BN-0. 

DBN-0 . 

RETURN 
END  IF 
END  IF 

COMPUTE  THE  BERNOULLI  FUNCTION  IN  THE  IirTERMEDlATE  RANGE 

IF (AX. CT. XLIM)  THEN 
ARG-DEXP (X) 


73 


noonnooooooooooononooo 


AKGl-AHG-1. 

BP-X/ARGl 

DBP-  (ABO*  (1  »*X)  -1 . )  /XMl/ARGl 

ARG-1 ./ABG 

ARGl-ABG-1. 

BN— X/XRGl 

DBN- (XBG* (1 .♦X) -1 . ) /ARCl/ABGl  ' 

ELSE 

C 

C  EVALUATE  THE  BERNOULLI  FUNCTION  FOR  SHALL  VALUES  OF  THE  ARGUMENT 


FP-l, 

FN-l. 

OF-1. 

GP-0 , 5 
GN-0 . 5 
DG-0.5 
SIGN-1. 

C 

C  EVALUATE  THE  TEFMS  OF  THE  SERIES 

C 

10  I-I'H 

AI-DFLOAT(I) 

SIGN— SIGN 
PF-DF*X/AI 

DG-DG*X*AI/ ( (Al+1 . ) * (Al-1 . )  ) 

FP-FP^DF 

rN-rN^DF*SICN 

GP-GP-*-DG 

CN-GN*DC*SICN 

IF(DABS (DG) .GT.tPS)  GO  TO  10 

BP-l./FP 

DBP— CP/  (FP*FP) 

BH-1 .  /f  t* 

DBN— GN/  (FNTN) 

END  IF 
C 

RETURN 

END 


SUBROUTINE  dbBrl  (X,BP,BN) 

»$$$S$$$$$$$$$$$$$$$$S$$S$$$S$$$$S$$$5$$$$$$5$$$S$$$$S$ 

Author  :  G.BAccaranifUniv.  of  Bologna 
Varaion  :  1.0 

Cat*  :  October  23  ,  1984 

^$$$S$$S$$S$$S$$$$$$$$$$3$$$$$$$$$$$$$SS$$$$$S$S$$S$$S$ 


This  routine  avaluataa  the  Bernoulli  function 
for  the  arguxnentf  X  and  -X 


Input  variable 

X  -  Argument  of  the  Bernoulli  function 

Output  variables: 
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BP 

BN 


Bernoulli  function  of  X 
Bornoulli  function  of  •X 


DATA: 


XLIM 

EPS 

KD 

KDIGIT 


10. ♦♦(-ND) 

10.** (-KDIGIT) 

nun\b«r  of  adnitt^d  arronaous  digits 
numbar  of  significant  digits 


IMPtlClT  DOUBLE  PRECISION  (A-H,0-Z) 
DATA  XLIM, EPS/0. OIDO,  X.D-15/ 


Exacutabla  coda. 


AX-DABS (X) 

Compute  the  asymptotic  valuas  of  the  Barnoulli  function. 

IF(AX.GT.80  DO)  THEN 
1F(X.GT.0 .DO)  THEN 
BP-0. DO 
BN-X 
RETURN 

ELSE 

BP— X 
BN-O.DO 

return: 

END  IF 
END  IF 

Compute  the  Barnoulli  function  in  the  intannadiata  range. 

IF (AX. CT. XLIM)  THEN 
ARG  -DEXP(X) 

ARGl-ARG-l.DO 
BP  -X/ARGl 
ARG  -1. DO /ARG 
ARGl-ARG-l.DO 
BN  —X/ARGl 

ELSE 


Evaluate  the  Barnoulli  function  for  small  valuas  of  tba  argumant , 
I-l 

FP-l.DO 

FN-l.DO 

DF-l.DO 

CP-O . 5DC 

GN-0.5DG 

DG-0.5D0 

SIGN  .DO 

Evaluate  the  terms  of  the  series. 


Al-DFLOAT(I) 
SIGN— SIGN 
DF-DF*X/AI 
FP-FP^DF 
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FN-FN+or-SIGN 
IF(DXBS(DF) .GT 
BP«1.00/FP 
BN-l.DO/FN 
END  XF 
RTTORN 
END 


.EPS)  CO  TO  10 
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B.3  dgtsl.f 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


•ubrcutln«  dgt si (n, c, d, a, b, inf o ) 
inta9ar  n^info 

doubla  precision  c(301) ,d(301) ,a(301) ,b(301) 

dOTsI  given  a  general  tridiagonal  oiatrix  and  a  right  hand 
aide  will  find  the  solution. 

on  entry 

n  integer 

is  the  order  of  the  tridiegonal  matrix. 

c  double  precision(n) 

is  the  subdaagonal  of  the  tridiagonal  matrix. 

Z's2]  through  c(r.}  should  contain  the  subdiagcnal. 
on  output  c  is  destroyed. 

d  double  precision (n) 

is  the  diagonal  of  the  tridiagonal  matrix, 

on  output  d  is  destroyed. 

e  double  precis ’.on  (n) 

is  the  superdiagonal  of  the  tridiagonal  matrix. 
e(l)  through  e(n-l)  should  contain  the  superdiagonai . 
on  output  a  is  destroyed. 

b  double  precision (n) 

is  the  right  hand  side  vector. 


on  return 

b  is  the  solution  vector. 

info  integer 

•  0  normal  value. 

«  k  If  the  k*th  element  of  the  diagonal  becomes 
exactly  zero,  the  subroutine  returns  when 
this  is  detected. 

Unpack,  this  version  dated  Ofl/14/78  . 
jack  dongarra,  argonne  national  laboratory. 

no  externals 
fortran  dabs 

internal  variables 

integer  k, kb, kpl, nml, nm2 
double  precision  t 

begin  block  permitting  ...exits  to  100 

info  -  0 
c(l)  -  d(l> 
nml  •  n  -  I 

if  (nml  .It.  X)  go  to  40 
d{l)  -  e(l) 
e(l)  «  O.OdO 
e{n)  -  O.OdO 


do  30  k  *  1,  runl 
kpl  -  k  +  1 

find  the  largest  of  the  two  rows 
1  f  (djibs  ( c*  ( ) )  .  It  .  ( r  ( ^  ) )  )  go  to  10 


O  Q 


interchange  row 

t  •  c(kpl) 
c()cpl)  -  c(k) 
c(k)  •  t 
t  -  d(kpl) 
d(kpl)  -  d(k) 
d(k)  -  t 
t  -  •(kpl) 

•  (kpl)  -  a(k) 

•  (k)  -  t 
t  -  b(kpl) 
b(kpl)  •  b(k) 
b(k)  •  t 

10  continue 

c 

c  zero  alamants 

c 

if  (c(k)  .n«.  O.OdO)  go  to  20 
info  -  k 

c  . axit 

go  to  100 

20  continue 

t  -  -c(kpl)/c(k) 
c(kpl)  •  d(kpl)  t*d(k) 
d(kpl)  •  aikpl)  ♦  t*a(k) 

•(kpl)  •  O.OdO 

b(kpl)  •  b(kpl)  ♦  t*b(k) 

30  continue 

40  continue 

if  (c(n)  .n«.  O.OdO)  go  to  50 
info  -  n 
go  to  90 
SO  continue 

c 

c  back  solve 

c 

nin2  -  n  -  2 
b  (n)  -  b  (n)  /c  (n) 
if  (n  .eq.  1)  go  to  80 

b(nml)  »  (b(nml)  -  d (nml) *b (n) ) /c (nml) 
if  (nn\2  .It.  1)  go  to  70 
do  60  kb  •  1,  nm2 
k  •  nm2  -  kb  1 

b(k)  -  (b(k)  -  d(k)*b'k+l;  -  e  (k)  *b  (k<^2) ) /c  (k) 
60  continue 

70  continue 

80  continue 

90  continue 
100  continue 
c 

return 

end 
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t 


Hatrixl *• 

procvot  dattt :  9/14/87 

process  t«mp«ratu«:o:  625  C 

process  pressure  .200  torr 

flows:  30/40/80  scon 

injttctoi  positions:  D/C/6. 875in  R  of  C 


iiafttr 

id 

20 

top 

center 

bottom 

i*ft 

right 

ave 

unit 

3908.00 

3919.00 

3889.00 

3945.00 

3864.00 

3905.00 

0.78 

26 

3911.00 

4000.00 

3982.00 

3927.00 

3990.00 

3962.00 

1.01 

35 

4046.00 

4003.00 

4005.00 

4048.00 

4012.00 

4022.80 

0.56 

45 

4062.00 

4013.00 

4022.00 

4053.00 

4027.00 

4035.40 

0.52 

55 

4134.00 

4079.00 

4118.00 

4127.00 

4125.00 

4116.60 

0.53 

65 

4216.00 

4183.00 

4237.00 

4219.00 

4223.00 

4215.60 

0.47 

75 

4400.00 

42';9.00 

4480.00 

4411.00 

4400.00 

4394.00 

1.65 

85 

4628.00 

4634.00 

4693.00 

4  r, 53. 00 

4630.00 

4647.60 

0.59 

95 

4658.00 

4679.00 

4844.00 

4695.00 

4668.00 

4708.80 

1.63 

105 

4648.00 

4668.00 

4821.00 

4703.00 

4650.00 

4698.00 

1.54 

115 

4488.00 

4478,00 

4646.00 

4621.00 

4471.00 

4540.90 

1.80 

124 

4205,00 

4217.00 

4271.00 

4257.00 

4177.00 

4225,40 

0.91 

130 

3926.00 

3929.00 

3961.00 

3952.00 

3908.00 

3935.20 

0.54 

•ve  unif  :  0.97 

center  ave:  4237.000  *♦-/-  6.76 

ave  ave  :  4262.092  +/-  7.07 


deposition  tiine: 

75.000  min. 

average 

wafer 

yrowth  rate 

id 

A/min 

20 

52.07 

26 

52.83 

35 

53.64 

45 

53.81 

55 

54.89 

65 

56.21 

75 

58,59 

85 

61.97 

95 

62.78 

105 

62.64 

115 

60.54 

124 

56.34 

130 

52.47 

80 


matrix2 .€ 

procftsa  date:  9/16/87 

process  temperature:  €25  C 

process  pressure  :  .200  torr 

silane  flows:  45/55/50  seem 

injector  positions:  D/C/80875in  R  of  C 


wafer 

id 

20 

top 

center 

bottom 

left 

right 

ave 

unif 

4C23.G0 

4019.00 

3984.00 

4053.00 

3950.00 

4005.80 

0.99 

26 

3966.00 

4051.00 

4023.00 

3975.00 

4028.00 

4008.60 

0.91 

35 

4040.00 

3996.00 

3990.00 

4034.00 

4002.00 

4012.40 

0.57 

45 

3978.00 

3928.00 

3936.00 

3968.00 

3951.00 

3952.20 

0.53 

55 

3928.00 

3898.00 

3921.00 

3914.00 

3924.00 

3917.00 

0.30 

€5 

3920.00 

3907.00 

3969.00 

3913.00 

3910.00 

3929.80 

0.64 

75 

3934.00 

3934.00 

4037.00 

3934.00 

3967.00 

3961.20 

1.13 

85 

3952.00 

3937.00 

4008.00 

3967.00 

3971,00 

3967.00 

0.67 

95 

3966.00 

3961.00 

4036.00 

3986.00 

3986.00 

3987.00 

0.74 

105 

3969.00 

3998.00 

4226.00 

4018.00 

3993.00 

4040.80 

2.60 

115 

3907.00 

3929.00 

4087.00 

.3950.00 

3933.00 

3961.20 

1.82 

124 

3819.00 

3835.00 

3886.00 

3851.00 

381S  00 

3842.00 

0.73 

130 

3389.00 

3396.00 

3437.00 

3413.00 

3371 .00 

3401.20 

0  .74 

sve  unif  0.95 

center  ave:  3906.846  +/-  4.19 

sve  ave  3922.015  ♦/•  4.19 


deposition  tim: 

75.000  min 

average 

wafer 

growth  rate 

id 

A/min 

20 

53.41 

26 

53.45 

35 

53.50 

45 

52.70 

55 

52.23 

€5 

52.40 

75 

52.82 

85 

52.89 

95 

53.16 

105 

53.88 

115 

52.82 

124 

51.23 

130 

45.35 

81 


matrix3 


process  date: 
process  temperature: 
process  pressure 
silane  flows;  60/70/20 
injector  positions: 


9/16/87 
625  C 
.200  torr 
ecczn 

D/C/10. 875tn  R  of  C 


wafer 

id 

top 

center 

bottom 

20 

4386.00 

4277.00 

4248.00 

26 

4206.00 

4279,00 

4264.00 

35 

4264.00 

4231.00 

4223.00 

45 

4159.00 

4096.00 

4104.00 

55 

4004.00 

3970,00 

4011.00 

65 

3939.00 

3927.00 

4008.00 

75 

3889.00 

3899.00 

4014.00 

65 

3826.00 

3817.00 

3844.00 

95 

3621.00 

3470.00 

3489.00 

105 

3397.00 

3352.00 

3427.00 

115 

3269.00 

3256.00 

3385.00 

124 

3128.00 

3138,00 

3230.00 

130 

3010.00 

3013.00 

3029.00 

i«ft 

right 

ave 

unif 

4288.00 

4207,00 

42£  20 

1.55 

4222.00 

4272.00 

4246.60 

0.77 

4263.00 

4232.00 

424.'  .60 

0.46 

4160.00 

4126.00 

4129.00 

0.72 

3997.00 

4011.00 

3998.60 

0.43 

3937.00 

3970.00 

3956.20 

0.84 

3885.00 

3924.00 

3922.20 

1.36 

3828.00 

3829.00 

3828.80 

0.25 

3621.00 

3609.00 

3562.00 

2.13 

3426.00 

3375.00 

3395.40 

0.96 

3326.00 

3254.00 

3298.40 

1.71 

3199.00 

3109, OC 

3160.80 

1.62 

3022.00 

2862.00 

2987.20 

2.36 

ave  unif  :  1.17 

center  ave:  3748.231  ♦/-  11.96 

ave  ave  :  3770.077  +/-  11.73 


deposition  time:  75.000  min 


wafer 

average 
growth  rate 

id 

A/min 

20 

57.08 

26 

56.65 

35 

56.57 

45 

.^5.05 

55 

53.31 

65 

52.75 

75 

52.30 

85 

51.05 

95 

47.49 

105 

45.27 

115 

43.98 

124 

42.14 

130 

39.83 

82 


4 


matrix4 . * 

process  date:  9/18/87 

process  temperature:  625  C 

process  pressure  :  .250  torr 

silane  flows:  30/55/65  seem 

injector  positions:  D/C/10. 875in  R  of  C 


wafer 

id 

top 

center 

bottom 

20 

3480.00 

3470.00 

3453.00 

26 

3453.00 

3488.00 

3487.00 

35 

3648.00 

3478.00 

3577,00 

45 

3635.00 

3465.00 

3475.00 

55 

3673.00 

3625.00 

3685.00 

€5 

3831.00 

3824.00 

3874.00 

75 

3886.00 

3883.00 

4011.00 

85 

3981.00 

3946.00 

4021.00 

95 

4058.00 

4026.00 

4112.00 

105 

4211.00 

4217.00 

4456.00 

115 

4241.00 

4271.00 

4672.00 

124 

4143.00 

4222.00 

4648,00 

130 

3948.00 

3983.00 

4100.00 

ave 

unif  : 

1.95 

center  ave: 

3838.308 

+/- 

7.96 

ave 

ave  : 

3886.538 

8.02 

left 

right 

ave 

unif 

3620. OC 

3432.00 

3401.00 

2.13 

3473.00 

3602.00 

3500.60 

1.67 

3642.00 

3593.00 

3587.60 

1.91 

3631.00 

3464.00 

3538.00 

2.46 

3675.00 

3687.00 

3669.00 

0.69 

3837.00 

3845.00 

3642.20 

0.50 

3889.00 

3906.00 

3915,00 

1.39 

3999.00 

3990.00 

3987.40 

0,69 

4088.00 

4056.00 

4068. OC 

0.81 

4272.00 

4209.00 

4273.00 

2.47 

4425.00 

4258.00 

4373,40 

4.17 

4260.00 

4148.00 

4284.20 

4.89 

4002.00 

3945.00 

3995.60 

1.58 

deposition  time:  75.000  min 


wafer 

average 
growth  rate 

id 

A/min 

20 

46.55 

26 

46.67 

35 

47.83 

45 

47.:: 

55 

48.92 

65 

51.:'' 

75 

52.20 

85 

53.17 

95 

54.24 

105 

56.97 

115 

58.31 

124 

57.12 

130 

53.27 

33 


matrixS .€ 

process  date:  9/18/87 

process  temperature:  625  C 

process  pressure  :  *250  torr 

silane  flows:  45/70/35  scon 

injector  positions:  D/C/6. 875in  R  of  C 


wafer 

id 

20 

top 

center 

bottom 

left 

right 

eve 

unif 

414'7.00 

4134.00 

4087.00 

4179.00 

4042.00 

4117.80 

1.31 

26 

4060.00 

4171,00 

4129.00 

4089.00 

4161.00 

4122.00 

1.14 

35 

4159.00 

4087.00 

4084.00 

4156.00 

4096.00 

4116.40 

0.92 

45 

4077.00 

4017.00 

4032.00 

4065.00 

4051.00 

4048.40 

0.60 

55 

4016.00 

3981.00 

4037.00 

4008.00 

4034.00 

4015.20 

0.56 

65 

4007.00 

4000.00 

4120.00 

4013.00 

4057,00 

4039,40 

1.24 

75 

4016.00 

4035.00 

4222.00 

4022.00 

4089.00 

4076.80 

2.11 

85 

3974.00 

3971.00 

4088.03 

4000.00 

4012.00 

4009.00 

1.18 

95 

3921.00 

3929.00 

4067,00 

3957.00 

3954.00 

3965.60 

1.48 

105 

3836.00 

3839,00 

3946.00 

3858.00 

3835.00 

3862.80 

1.23 

115 

3482.00 

3464.00 

3648.00 

3620.00 

3481.00 

3539.00 

2.47 

124 

3267.00 

3260,00 

3327.00 

3315.00 

3256.00 

3265.00 

1.02 

130 

3052.00 

3043.00 

3067.00 

3056.00 

3040,00 

3051,60 

0.35 

eve  unif 

1.20 

center  ave: 

3840.846 

+/- 

9.24 

ave  ave 

3865.308 

+  /- 

9.02 

deposition  time: 

75.000  min 
average 

wafer 

growth  rate 

id 

A/min 

20 

54.90 

26 

54.96 

35 

54.89 

45 

53.98 

55 

53.54 

65 

53.86 

75 

54.. ’6 

85 

53.45 

95 

52.87 

105 

51.50 

115 

47. la 

124 

43.80 

130 

40.69 

84 


matrixS .e 

process  date:  9/21/87 

process  temperature:  625  C 

process  pressure  :  .250  torr 

silane  flows:  60/40/50  scan 

injector  positions:  D/C/8, 875in  R  of  C 


wafer 

id 

top 

center 

bottom 

20 

4550.00 

4547.00 

4494.00 

26 

4400.00 

4529.00 

4497.00 

35 

4486.00 

4404.00 

4393.00 

45 

4307.00 

4213.00 

4215.00 

55 

4106.00 

4051.00 

4085.00 

65 

4056.00 

4017.00 

4080.00 

75 

4044.00 

4025.00 

4128.00 

85 

4037.00 

4007.00 

4074.00 

95 

4043.00 

4025.00 

4125.00 

105 

4032.00 

4073.00 

4418.00 

115 

3953.00 

3989.00 

4201.00 

124 

3819.00 

3844.00 

3922.00 

130 

3420.00 

3427.00 

3479.00 

ave 

unif  : 

1.23 

center  ave: 

4088.538 

7.23 

ave 

ave  : 

4110.354 

+/- 

6.89 

left 

right 

ave 

unif 

4574.00 

4444.00 

4521.80 

1.16 

4428.00 

4500.00 

4470.80 

1.21 

4470.00 

4424.00 

4435.40 

0.92 

4268.00 

4249.00 

4250.40 

0.92 

4091.00 

4104.00 

4087.40 

0.54 

4046,00 

4063.00 

4052.40 

0.58 

4041.00 

4062.00 

4060.00 

0.99 

4044.00 

4042.00 

4040.80 

0.59 

4050.00 

4054.00 

4061.00 

0.94 

4079.00 

4064.00 

4133.20 

3.88 

4005.00 

3988.00 

4027.20 

2.46 

3860,00 

3830.00 

3855.00 

1.05 

3445.00 

3425.00 

3439.20 

0.70 

deposition  time:  75.000  min 


wafer 

average 
growth  rate 

id 

A/min 

20 

60.29 

26 

59.61 

35 

59.14 

45 

56.67 

55 

54.50 

65 

54.03 

75 

54.13 

85 

53. B8 

95 

54.15 

105 

55.11 

115 

53.70 

124 

51.40 

130 

45.86 

85 


matrix^ .e 

process  date:  9/21/87 

process  temperature:  €25  C 

process  pressure  .350  torr 

silane  flows:  30/70/50  seem 

injector  positions:  D/C/8. 875in  R  of  C 


wafer 

id 

20 

top 

center 

bottom 

l*£t 

right 

ave 

unif 

3686.00 

3661.00 

3631.00 

8710.00 

3614.00 

3660.40 

1.07 

26 

3624.00 

3693.00 

3684,00 

3655.00 

3713.00 

3673.80 

0.95 

35 

3822.00 

3667,00 

3681.00 

3820.00 

3712.00 

3740.40 

2.01 

45 

3736.00 

3659.00 

3701.00 

3733.00 

3721.00 

3710.00 

0.85 

55 

3846.00 

3836.00 

3883.00 

3854.00 

3875.00 

3858.80 

0.51 

65 

3901.00 

3903,00 

4042. CO 

3906.00 

3982.00 

3946.80 

1.60 

75 

3996.00 

4024.00 

4274.00 

4016.00 

4095.00 

4081.00 

2.80 

85 

4059.00 

4057,00 

4216.00 

4095.00 

4129.00 

4111.20 

1.60 

95 

4092.00 

4090.00 

4242.00 

4136.00 

4153.00 

4142.60 

1.50 

105 

4103.00 

4170.00 

4655.00 

4210.00 

4190.00 

4265.60 

5.19 

115 

4010.00 

4063.00 

4437.00 

4087.00 

4068.00 

4133.00 

4.17 

124 

3864.00 

3883.00 

4015.00 

3906.00 

3877.00 

3909.00 

1.57 

130 

3449.00 

3449.00 

3484.00 

3471.00 

3461.00 

3462.80 

0.43 

ave  unif  ;  1.86 

center  ave:  3858.077  +/-  5,66 

ave  ave  3899.646  ♦/-  6.15 


deposition  time: 

75.000  min 

average 

wafer 

growth  rate 

id 

A/min 

20 

48.81 

26 

48.98 

35 

49.87 

45 

49.47 

55 

51.45 

65 

52.62 

75 

54.41 

85 

54.82 

95 

55.23 

105 

56.87 

115 

55.11 

124 

52.12 

130 

46.17 

matrixO  .e 

process  date:  9/23/87 

process  temperature:  625  C 

process  pressure  :  .350  torr 

silane  flows:  45/40/65  seem 

injector  positions:  D/C/10. 875in  R  of  C 


wafer 

id 

20 

top 

center 

bottom 

left 

right 

ave 

unif 

4105.00 

4093.00 

4057.00 

4135.00 

4031.00 

4084.20 

1.00 

26 

3978.00 

^.051.00 

4031.00 

4000.00 

4044.00 

<C 10.80 

0.77 

35 

4013.00 

3955.00 

3955.00 

4013.00 

3976.00 

5932.40 

0.73 

45 

3921.00 

3876.00 

3881.00 

3918.00 

3901.00 

3C99.40 

0.53 

55 

3872.00 

3845.00 

3869.00 

3872.00 

3873.00 

3866.20 

0.31 

65 

3876.00 

3861.00 

3925.00 

3896.00 

3896.00 

3890.80 

0.62 

75 

3933.00 

3916,00 

4042.00 

3940.00 

3961.00 

3953.40 

1.25 

85 

4032,00 

3996.00 

4075.00 

4050.00 

4036.00 

4037.80 

0.71 

95 

4132.00 

4097.00 

4195.00 

4184.00 

4127.00 

4147.00 

0.99 

105 

4286.00 

4408.00 

4414.00 

4492.00 

4396.00 

4399.20 

1.68 

115 

4440.00 

4624.00 

4855.00 

4641.00 

4476.00 

4607.20 

3.57 

124 

4260.00 

4462.00 

4784.00 

4620.00 

4282.00 

4481.60 

4.99 

130 

4045.00 

4096.00 

4247.00 

4134.ee 

4064.00 

4117.20 

1.94 

ave  unif  :  1.47 

center  ave:  4098.462  +/-  6.06 

ave  ave  ;  4114.785  +/-  5,76 


deposition  time: 

75.000  min 

average 

wafer 

growth  rate 

id 

A/min 

20 

54.46 

26 

55.61 

35 

53.10 

45 

51.99 

55 

51.55 

65 

51.86 

75 

52.78 

85 

53.84 

95 

55.29 

105 

58.66 

115 

61.43 

124 

59.75 

130 

54.90 

87 


inatrix9 

.e 

process 

date: 

9/23/87 

process 

temperature : 

623  C 

process 

pressure 

: 

.350  torr 

silene 

flows : 

60/55/35 

scan 

injector  positions: 

D/C/6. 87Sin 

R  Of  C 

wafer 

4  A 

top 

center 

bottom 

left 

right 

eve 

unlf 

io 

20 

4742.00 

4734.00 

4701.00 

4768.00 

4692.00 

4727.40 

0.66 

26 

4665.00 

4699.00 

4684.00 

4676.00 

4700.00 

4684.80 

0.32 

35 

4670.00 

4637.00 

4637.00 

4667.00 

4651.00 

4652.40 

0.34 

45 

4273.00 

4238.00 

4246.00 

4273.00 

4267,00 

4259.40 

0.38 

55 

4174.00 

4097.00 

4170.00 

4167.00 

4181.00 

4157.80 

0.83 

65 

4109,00 

4060.00 

4189.00 

4121.00 

4158.00 

4127.40 

1.19 

75 

4063.00 

4059.00 

4239.00 

4088.00 

4160.00 

4121.80 

1.87 

85 

3989.00 

3967.00 

4121.00 

4019.00 

4012.00 

4021.60 

1.47 

95 

3926.00 

3927.00 

4124.00 

3988.00 

3942.00 

3981.40 

2.10 

105 

3842.00 

3843.00 

3968.00 

3873.00 

3850.00 

3875.20 

1.38 

115 

3415.00 

3406.00 

3467.00 

3452.00 

3417.00 

3431.40 

0.77 

124 

3109.00 

3170.00 

3246.00 

3232.00 

3181.00 

3203.60 

1.04 

130 

3033.00 

2668.00 

3037.00 

3036.00 

n3032.00 

2961.20 

5.54 

ave 

unif  : 

1.38 

center  ave: 

3961.923 

♦/- 

15.17 

ave 

ave 

4015.800 

+/- 

13.63 

deposition  time: 

75.000  min 

average 

wafer 

growth  rate 

id 

A.^min 

20 

63.03 

26 

62.46 

35 

62.03 

45 

56.79 

55 

55.44 

65 

55.03 

75 

54.96 

85 

53.62 

95 

53.09 

105 

51.67 

115 

45.75 

124 

42.71 

130 

39.48 

88 


optixTiizeci  parameter  run 
process  date :  X2/A/B1 

process  temperature:  625  C 

process  pressure  :  .200  torr 

silane  flows:  45/40/65  seem 

injector  positions:  D/C/8. 875in  R  of  C 


wafer 

id 

20 

top 

center 

bottom 

left 

right 

ave 

unif 

3905.00 

3887.00 

3851.00 

3976.00 

3824,00 

3888.60 

1.49 

26 

3872.00 

3874 .00 

3862. GO 

3820.00 

3930. CO 

3873.20 

0.95 

35 

3885.00 

3844.00 

3849.00 

3895.00 

3855.00 

3665.60 

0.59 

45 

3861.00 

3826.00 

3833.00 

3856.00 

3845.00 

3844.20 

0.39 

55 

3835.00 

3700.00 

3830.00 

3842.00 

3837,00 

3808.80 

1.60 

€5 

3852.00 

3827.00 

3872.00 

3852.00 

3856.00 

3851.80 

0.42 

75 

3083.00 

3869.00 

3940.00 

3884.00 

3895.00 

3894.20 

0,70 

85 

3920.00 

3899.00 

3968.00 

3962.00 

3920.00 

3933.80 

0.76 

95 

3973.00 

3960,00 

4086.00 

4035.00 

3960.00 

4002.80 

1.40 

105 

3956.00 

4010.00 

4270.00 

4065.00 

3968,00 

4053.80 

3.16 

115 

3889.00 

3804.00 

4041.00 

3985.00 

3892,00 

3940.20 

1.76 

124 

3879.00 

3822.00 

3854.00 

3850.00 

3880.00 

3857.00 

0.62 

130 

3325.00 

3312.00 

3366.00 

3400.00 

3285.00 

3337.60 

1.36 

ave  unif  :  1.17 

center  ave:  3824.923  +/-  4.47 

ave  ave  :  3857.815  +/-  4.41 


deposition  time: 

75.000  min 

average 

wafer 

growth  rate 

id 

A/min 

20 

51.85 

26 

51.64 

35 

51.54 

45 

51.26 

55 

50.78 

65 

51.36 

75 

51.92 

85 

52.45 

95 

53.37 

105 

54.05 

115 

52.54 

124 

51.43 

130 

44.50 

89 


